REPORT  DOCUMENTATION  PAGE 

_ _ _ _ _ J _ _ 1 _ 

Public  reporting  burden  for  this  collection  of  infomnation  is  estimated  to  average  1  hour  per  response,  including  the  tim 
and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  Information.  Send  comitents  reg 
Information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  infoi 
1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  management  and  Budget,  Paperwor1<  Reduction  Project  (0704-018 

1.  AGENCY  USE  ONLY  (Leave  Blank)  1 2.  REPORT  DATE  1 3.  REPORT  TY 


Form  Approved 

OAvfO  Ai-  ‘  ' 


;^FRL-SR-BL-TR-9*- 


4.  TITLE  AND  SUBTITLE 

Application  of  a  Genetic  Algorithm  to  the  Optimization  of  a  Missile  Autopilot 
Controller  for  Performance  Criteria  with  Non-Analytic  Solutions _ 


6.  AUTHORS 

Richard  A.  Hull 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Central  Florida 


5.  FUNDING  NUMBERS 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AFOSR/NI 

110  Duncan  Avenue,  Room  B-115 

Bolling  Air  Force  Base,  DC  20332-8080  _ 


11.  SUPPLEMENTARY  NOTES 


10.  SPONSORING/MONITORING 
,  AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Approved  for  Public  Release 


13.  ABSTRACT  (Maximum  200  words) 

See  attached. 


12b.  DISTRIBUTION  CODE 


"Bestir  1C 

Approved  for  puLiic 
Distribution  Uniiniiied 


19980505  050 


15.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

Unclassified  Unclassified  Unclassified  UL 


Standard  Form  298  (Rev,  2-89) 

Prescribed  by  ANSI  Std.  239.18 
Designed  using  WordPerfect  6.1.  AFOSR/XPP.  Oct  96 


Abstract 


Modern  optimal  control  theory  provides  analytic  solutions  for  a  set  of  lin¬ 
ear  feedback  design  problems  with  linear  quadratic  performance  criteria.  Recent 
progress  in  the  field  of  robust  multivariable  feedback  design  has  incorporated  addi¬ 
tional  constraints  which  have  addressed  the  classical  concerns  with  stability  margins, 
system  sensitivity  and  disturbance  rejection.  Despite  these  important  advances, 
many  practical  design  problems  arise  in  which  the  desired  system  performance  con¬ 
straints  cannot  be  accommodated  by  the  available  theoretic  techniques. 

Genetic  algorithms  {  GA’s  ),  on  the  other  hand,  offer  a  numerical  search 
method  which  does  not  require  a  statement  of  the  mathematical  relationship  be¬ 
tween  the  performance  criteria  and  the  parameter  update  rule.  The  objective  of  this 
thesis  is  to  demonstrate  that  GA’s  provide  a  method  of  optimizing  control  system 
problems  with  analytically  intractable  constraints. 

A  linear  missile  airframe  and  actuator  state  space  model  is  developed  with 
linear  feedback  controller,  and  implemented  in  a  discrete  time  simulation.  A  genetic 
algorithm  is  constructed  to  optimize  the  linezu-  controller  parameters,  first  with 
respect  to  a  weighted  linear  quadratic  performace  index.  Additional  performance 
constraints  are  then  imposed  to  meet  rise  time,  peak  actuator  effort,  and  settling 
error  specifications.  Computer  simulation  results  show  that  the  genetic  algorithm 
provides  good  convergence  to  near  optimal  controller  designs  for  each  successive 
combination  of  constraints. 


Application  of  a  Genetic  Algorithm  to  the  Optimization  of  a  Missile 
Autopilot  Controller  for  Performance  Criteria  with  Non-Analytic 

Solutions 


by 

Richard  A.  Hull 

B.S.  E.S.M.,  University  of  Florida,  1972 


THESIS 

Submitted  in  partial  fulfillment  of  the  requirements 
for  the  degree  of 

Master  of  Science  in  Electrical  Engineering 
College  of  Engineering 
University  of  Central  Florida 
Orlando,  Florida 


Spring  Term 
1993 


r' 


' '  I 

,  ..i- 


;  ?  f  ‘  i  ’ 


/  r  ^ 

i  ■>'-  M.  I 


r.:>?  ;-,-wiewGd  and  is 

Iaw  AFR  190-12 

■^;L  <o;'!  in  ttd. 

Approvsd  f  oT  pnb''.  5  o  r»'f  - 
*3i2trib-.:tion  UiiliK-;-;  cd. 


'>  ■. 
--  *• 


0  l^rograra  Manage^ 


Abstract 

Modern  optimal  control  theory  provides  analytic  solutions  for  a  set  of  lin¬ 
ear  feedback  design  problems  with  linear  quadratic  performance  criteria.  Recent 
progress  in  the  field  of  robust  multivariable  feedback  design  has  incorporated  addi¬ 
tional  constraints  which  have  addressed  the  classical  concerns  with  stability  margins, 
system  sensitivity  and  disturbance  rejection.  Despite  these  important  advances, 
many  practical  design  problems  arise  in  which  the  desired  system  performance  con¬ 
straints  cannot  be  accommodated  by  the  available  theoretic  techniques. 

Genetic  algorithms  (  GA’s  ),  on  the  other  hand,  offer  a  niunerical  search 
method  which  does  not  require  a  statement  of  the  mathematical  relationship  be¬ 
tween  the  performance  criteria  and  the  parameter  update  rule.  The  objective  of  this 
thesis  is  to  demonstrate  that  GA’s  provide  a  method  of  optimizing  control  system 
problems  with  analytically  intractable  constraints. 

A  linear  missile  airframe  and  actuator  state  space  model  is  developed  with 
linear  feedback  controller,  and  implemented  in  a  discrete  time  simulation.  A  genetic 
algorithm  is  constructed  to  optimize  the  linear  controller  parameters,  first  with 
respect  to  a  weighted  linear  quadratic  performace  index.  Additional  performance 
constraints  are  then  imposed  to  meet  rise  time,  peak  actuator  effort,  and  settling 
error  specifications.  Computer  simulation  results  show  that  the  genetic  algorithm 
provides  good  convergence  to  near  optimal  controller  designs  for  each  successive 
combination  of  constraints. 
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CHAPTER  1 


INTRODUCTION 


Classical  autopilot  design  methods  rely  on  frequency  domain  techniques  to  achieve 
gain  and  phase  margin  criteria,  while  providing  desired  characteristics  of  the  closed 
loop  system  response  such  as  minimum  rise  time  and  maximum  overshoot. 

Autopilot  design  methods  based  on  modern  state  space  concepts  have  been 
developed  which  allow  the  designer  to  shape  the  dynamic  response  of  the  closed 
loop  system  by  placing  the  poles  of  the  compensated  system,  provided  issues  of 
controllability  and  observabibility  have  been  adaquately  addressed.  Determining 
where  to  place  the  poles,  however,  is  somewhat  more  subjective. 

Recent  progress  in  the  field  of  robust  multivariable  feedback  design  theory 
has  employed  singular  value  decomposition  as  a  measure  of  system  performance 
and  robustness.  Singular  value  loop  shaping  has  then  been  accomplished  via  Hoo  , 
frequency  weighted  LQG,  and  LQG  loop  transfer  recovery  design  techniques  [11,4]. 

All  of  these  important  advances,  however,  have  placed  constraints  on  the 
ways  in  which  system  performance  and  robustness  characteristics  must  be  specified. 
Most  control  system  engineers  will  be  faced  with  many  design  problems  in  which  the 
given  specifications  cannot  readily  be  formulated  in  terms  required  by  the  available 
theoretic  design  criteria.  In  particular,  many  long  accepted  norms  of  classical  system 
performance  criteria,  such  as  rise  time  and  maximum  overshoot,  may  result  in 
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problem  formulations  which  do  not  have  analytic  solutions. 

Genetic  algorithms  (GA’s),  on  the  other  hand,  offer  a  numerical  search 
method  which  mathematically  mimics  some  of  the  aspects  of  biological  natural 
selection  in  order  to  achieve  adaptation  to  a  specific  measure  of  fitness.  An  impor¬ 
tant  advantage  of  genetic  algorithms  is  that  they  do  not  require  a  mathematical 
statement  of  the  relationship  between  the  performance  criteria  and  the  parameter 
update  rule.  Instead  they  rely  on  reproductive  characteristics  of  a  sample  gene  pool 
which  is  bred  to  maximize  the  desired  fitness  measure.  Furrthermore,  genetic  algo¬ 
rithms  are  good  at  finding  global  as  opposed  to  local  optimums.  These  properties 
make  GA’s  an  excellent  candidate  for  numerically  optimizing  a  multicriterion  design 
specification  problem  which  incorporates  an  arbitrary  but  feasible  combination  of 
performance  and  constraint  criteria. 

The  objective  of  this  thesis  is  to  demonstrate  that  a  genetic  algorithm  pro¬ 
vides  a  method  for  optimizing  an  autopilot  control  system  with  respect  to  a  set  of 
design  specifications  which  are  feasible  but  have  an  analytically  intractable  solution. 


1.1  The  Control  System  Optimization  Problem 

The  controller  design  problem  begins  with  a  model  of  the  system  to  be  controlled, 
referred  to  as  the  plant  ,  and  a  set  of  design  goals,  or  specifications  .  The  plant  in 
this  instance  will  be  regarded  as  the  physical  system  (i.e.  missile  airframe)  along 
with  its  control  surface  actuators  and  feedback  sensors.  The  design  specifications 
must  incorporate  criteria  related  to  the  performance,  stability,  and  robustness  of 
the  resulting  controlled  system. 

In  the  general  case,  the  form  of  the  controller,  as  well  as  the  form  of  the 
design  specifications,  may  all  be  open  for  selection  by  the  designer.  In  fact,  in 
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the  most  general  case,  decisions  regarding  the  design  of  the  plant,  the  number  and 
fidelity  of  feedback  sensors,  the  size  and  location  of  control  surfaces,  the  bandwidth 
of  actuators,  and  the  like,  may  also  be  determined  or  influenced  by  the  control 
systems  engineer.  At  some  stage  in  the  design  process,  however,  the  other  factors 
become  fixed,  and  the  controller  design  problem  becomes  a  search  for  a  suitable 
controller  that  will  satisfy  the  stated  design  specifications. 

In  optimal  control  theory,  the  goal  of  finding  the  best  controller  is  introduced. 
To  this  goal,  the  design  specifications  are  formulated  in  mathematical  terms  as 
a  performance  measure  plus  auxiliary  constraints.  The  objective  of  the  optimal 
control  problem  is  to  determine  the  controller  that  will  minimize  (or  maximize)  the 
performance  measure,  and  at  the  same  time  satisfy  the  auxiliary  constraints  [9]. 

In  practice,  most  optimal  control  procedures  fix  the  form  of  the  controller 
and  codify  the  form  of  the  performance  measure  in  order  to  develop  analytical  solu¬ 
tions.  Most  optimal  controls  work  is  based  on  Linear  Time  Invariant  Causal  (LTIC) 
Systems  with  linear  state  feedback  controllers.  The  best  known  of  these  methods 
employs  the  weighted  mean-square  error  of  the  regulation  variables  and  control 
effort  to  form  a  quadratic  performance  index.  This  usually  results  in  the  Linear 
Quadratic  Regulator  (LQR)  problem  when  full  state  feedback  is  available,  or  the 
Linear  Quadratic  Gaussian  (LQG)  problem  when  state  estimation  is  required.  The 
LQR  design,  while  extremely  robust  from  a  theoretical  standpoint,  is  often  imprac¬ 
tical  to  implement,  and  the  LQG  design  results  in  a  controller  that  is  notoriously 
non-robust. 

In  their  book  Linear  Controller  Design,  Limits  of  Performance  ,  Boyd 
and  Barratt  [2]  develop  an  analytical  framework  for  determining  the  most  stringent 
set  of  controller  design  specifications,  within  a  specific  class,  that  can  be  met  using 
any  controller  design  method  for  a  given  plant  and  feedback  control  configura- 
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tion.  This  results  in  the  so-called  Pareto  optimal  specification  ,  and  the  resulting 
controller  is  the  Pareto  optimal  design  .  After  classifying  many  types  of  design 
specifications,  Boyd  and  Barratt  present  a  numerical  method  for  solving  controller 
design  problems  with  closed-loop  convex  nondifferentiable  design  specifications  that 
do  not  have  analytic  solutions.  While  this  is  an  impressive  and  inspiring  body  of 
work,  such  Pareto  optimal  designs  usually  result  in  controllers  of  inordinatly  high 
order  (say  100  or  more)  for  even  relatively  simple  systems  ! 

For  the  purpose  of  this  thesis,  the  control  system  optimization  problem  will  be 
limited  to  a  fixed  form  of  controller:  linear  state  feedback  control;  and  a  given  Linear 
Time  Invariant  Causal  (LTIC)  plant.  A  performace  measure  and  several  auxiliary 
constraints  will  be  successively  incorporated  into  a  single  fitness  measure  that  can 
be  maximized  by  a  genetic  algorithm.  It  will  be  demonstrated  that  the  genetic 
algorithm  can  be  successfully  applied  to  the  resulting  multicriterion  optimization 
problems  for  which  analytic  solutions  do  not  exist. 

1.2  Common  Measures  of  Performance 

Control  system  design  specificatons  must  incorporate  criteria  related  to  the  per¬ 
formance,  stability,  and  robustness  of  the  resulting  controlled  system.  In  classical 
controls  design  methods,  performance  specifications  describe  how  the  closed  loop 
system  should  respond  to  commanded  inputs  and  how  well  it  should  reject  un¬ 
commanded  disturbances.  Some  examples  of  classical  performance  specifications 
include:  rise  time,  percent  overshoot,  and  settling  time  [14]. 

Robustness  specifications  are  intended  to  limit  the  variations  in  performance 
of  the  closed  loop  system  due  to  variations  in  the  plant  during  operation,  or  due 
to  errors  and  simplifications  in  the  plant  model  used  for  design.  In  classical  de- 
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sign  methods,  robustness  specifications  usually  take  the  form  of  allowable  gain  and 
phase  margins.  Closed  loop  system  stability  is  usually  understood  to  be  a  hard 
constraint  that  is  contained  within  the  performance  specifications  and  guaranteed 
by  the  robustness  specifications. 

Modern  optimal  control  theory  provides  an  analytical  solution  to  the  Linear 
Quadratic  Regulator  (LQR)  problem  by  minimizing  a  quadratic  performance  index 
of  the  form  : 


J=f  [x'(t)Qx{t) +u'(t)Ru(t)]dt  (1.1) 

JO 

where  : 

x(t)  is  the  state  variable  vector, 
u(t)  is  the  control  variable  vector, 

Q  is  a  positive  semi-definite  weighting  matrix 
R  is  a  positive  definite  weighting  matrix 
T  is  the  finite  final  time 

For  the  dynamic  process  characterized  by  : 

x(i)  =  Ax(i)  +  Bu(t)  (1.2) 

solution  of  the  Algebraic  Riccati  Equation  (ARE)  will  provide  a  steady  state  linear 
control  law  of  the  form  : 

u{t)  =  -kx(t)  (1.3) 

where  k  is  the  steady  state  gain  matrix  that  minimizes  the  quadratic  performance 
index  J.  This  basic  perfomance  index  can  be  modified  to  incorporate  niimerous 
variations  such  as: 


1.  final  value  of  state  vector  at  time  T 
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2.  time  varying  Q  matrix 

3.  time  varying  R  matrix 

4.  non-zero  initial  time 

5.  infinite  final  time  {T  =  oo) 

1.3  Standard  Optimization  Techniques 

The  solution  of  the  control  system  optimization  problem  depends  to  a  large  extent 
upon  certain  properties  of  the  design  specification.  In  a  few  cases  closed  form 
analytical  solutions  may  exist,  but  usually  a  form  of  numerical  search  technique  is 
required.  Many  forms  of  search  techniques  have  been  documented,  but  most  fall 
into  one  of  three  broad  categories:  calculus-based,  enumerative,  or  random  [7]. 

If  the  design  specification  is  differentiable  ,  then  several  standard  analytical 
as  well  as  calculus-based  search  techniques  may  apply.  If  the  design  specification 
is  not  differentiable,  but  is  closed-loop  convex  [2,  chapter  6]  then  numerical  search 
methods  such  as  the  Nelder-Meade  algorithm  [12],  the  cutting  plane  algorithm 
[2],  or  the  ellipsoid  algorithm  [2]  may  be  required.  If  the  design  specification  is 
non-differentiable  and  non-convex,  then  the  more  time  consuming  enumerative  or 
random  search  techniques  may  be  the  only  hope. 

The  quadratic  performance  index  shown  in  equation  1.1  is  a  differentiable 
specification,  and  the  Linear  Quadratic  Regulator  problem  has  been  analytically 
solved  in  its  many  forms  by  application  of  the  Hamilton- Jacobi  equation  as  devel¬ 
oped  in  either  [1,  chapter  2]  or  [9,  chapter  5,  section  2].  The  solution  of  the  resulting 
two-point  boundary  value  problem  results  in  an  optimal  control  of  the  form: 


u*{t)  =  -k{t)x{t)  =  -R-'(f)B'(t)P(t)x(t) 


(1.4) 
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where  P(t)  is  the  solution  to  the  familiar  Riccati  equation: 

T{t)  =  -P(t)A(t)  -  A'{t)P(t)  -  Q(t}  +  P(t)B(t)R-^(OB'(t)P(t)  (1.5) 

The  optimal  LQR  control  is  thus  time  varying,  and  the  Riccati  equation 
can  be  solved  analytically  for  simple  systems,  or  numerically  for  more  complicated 
systems  by  integrating  backward  from  the  final  time  T.  Methods  of  solving  the 
Riccati  equation  are  presented  in  many  standard  texts,  see  [9,1,5,6,3]  for  details. 

The  steady  state  solution  of  equation  1.5  can  be  determined  by  setting  the 
derivative  to  zero  and  solving  the  resulting  Algebraic  Riccati  Equation  (ARE): 

0  =  -PA  -  A'P  -  Q  +  PBR-^B'P  (1.6) 

It  has  been  shown  that  if  the  system  is  asymptotically  stable  and  is  both  control¬ 
lable  and  observable,  then  the  ARE  has  a  unique  positive  definite  solution  P,  which 
minimizes  the  performance  index  J,  when  used  in  the  optimum  control  law  of  equa¬ 
tion  1.4.  This  approach  is  often  used  to  generate  a  non-time  varying  control  law 
that  produces  near  optimal  results. 

The  LQR  design,  while  extremely  robust  from  a  theoritical  standpoint,  re¬ 
quires  full  state  feedback  which  is  often  impractical  in  the  implementation  or  actual 
systems.  Optimal  state  estimation,  the  dual  problem  to  the  optimal  control  formu¬ 
lation,  usually  results  in  the  implementation  of  a  Kalman  filter  to  provide  estimates 
of  missing  or  noisy  states  in  the  feedback  path.  This  configuration  is  known  as 
leading  to  the  Linear  Quadratic  Gaussian  (LQG)  controller.  In  practice  however, 
the  LQG  design  results  in  a  controller  that  is  notoriously  non-robust. 

The  Nelder-Meade  algorithm,  is  a  numerical  search  technique  that  does  not 
require  analytical  differentiation  [12].  It  is  effective  and  computationally  compact, 
however,  it  is  not  guaranteed  to  find  global  optimums.  Instead,  it  must  be  initial¬ 
ized  with  an  initial  guess,  and  tends  to  find  the  local  minimum  of  a  function,  in  the 
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vicinity  of  the  initial  guess.  The  Nelder-Meade  algorithm  is  formulated  in  terms  of 
a  simplex,  which  is  a  generalized  triangle  of  N  dimensions.  For  functions  of  two  vari¬ 
ables,  the  simplex  is  just  a  triangle.  The  algorithm  works  by  evaluating  the  function 
to  be  optimized  at  the  three  vertices  of  the  triangle,  ranking  the  vertices  from  best 
to  worst.  The  worst  vertex  is  rejected,  and  a  method  is  provided  for  computing 
a  new  vertex  so  that  a  sequence  of  triangles  is  formed  that  will  converge  on  the 
minimum  of  the  function.  The  simplex  concept  easily  extends  the  mechanization  of 
the  algorithm  to  include  functions  of  N  variables. 

Modern  control  theory  optimization  techniques,  are  limited  in  the  ways  in 
which  system  performance  and  robustness  characteristics  must  be  specified.  Most 
control  system  engineers  will  be  faced  with  many  design  problems  in  which  the 
given  specifications  cannot  readily  be  formulated  in  terms  required  by  the  avail¬ 
able  theoretic  design  criteria.  In  particular,  many  long  accepted  norms  of  classical 
system  performance  criteria,  such  as  maximum  overshoot  or  settling  time,  may  re¬ 
sult  in  relationships  which  do  not  have  analytic  solutions.  In  many  other  situations, 
desirable  performance  or  robustness  criteria  may  result  in  relationships  between  con¬ 
troller  parameters  and  the  performance  index  which  are  either  non-differentiable  or 
non-convex.  Either  situation  makes  solution  by  standard  numerical  methods  diffi¬ 
cult.  In  particular,  numerical  search  algorithms  tend  to  get  stuck  on  local  maxima 
or  minima,  while  enumerative  or  random  search  techniques  take  too  long. 

1.4  An  Introduction  to  Genetic  Algorithms 

Genetic  algorithms  are  a  numerical  search  method  which  mathematically  mimics 
some  of  the  aspects  of  biological  natural  selection  in  order  to  achieve  adaptation  to  a 
specific  measure  of  fitness.  Genetic  algorithms,  or  GA’s,  first  appeared  in  the  early 
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1970’s  and  were  largly  developed  by  John  Holland,  his  colleagues  and  students  at  the 
University  of  Michigan  [7].  Being  non-calculas  based,  GA’s  do  not  require  derivative 
information  about  the  design  specification.  Instead  they  use  a  biological  approach 
to  randomly  mating  individual  genetically  coded  characteristics,  with  preference 
weighted  toward  those  which  produce  individuals  displaying  the  greatest  fitness  . 
While  employing  random  choice  as  a  tool  in  the  process,  genetic  algorithms  are 
not  purely  random  search  techniques,  but  rather  exploit  the  surprisingly  powerful 
optimization  tendency  of  population  group  evolution.  Of  particular  importance  is 
the  fact  that  genetic  algorithms  work  well  in  finding  global  optimums  and  are  much 
less  likely  to  converge  on  local  maxima  or  minima  in  the  manner  of  many  numerical 
search  methods. 

Genetic  algorithms  may  be  constructed  of  varying  degrees  of  complexity,  but 
an  effective  form  can  be  fairly  easy  to  implement.  First,  a  fitness  function  to  be 
maximized  by  the  algorithm  must  be  determined  for  the  problem  at  hand.  In  the 
case  of  an  LQR  control  system,  for  example,  this  fitness  function  could  be  related 
to  the  quadratic  performance  index.  In  this  case  the  fitness  function  must  express 
an  inverse  relationship  to  the  quadratic  performance  index,  since  the  LQR  solution 
works  to  minimize  this  index,  while  the  genetic  algorithm,  for  technical  reasons  that 
will  become  apparent,  works  to  maximize  the  fitness  function. 

Once  the  fitness  function  is  determined,  the  variable  parameters,  such  as 
the  control  system  feedback  gains,  are  coded  into  bit  strings  consisting  of  ones  and 
zeros.  An  initial  population  of  random  strings  is  generated,  and  the  system  fitness 
function  is  evaluated  for  each  individual  in  the  population.  Successive  generations 
are  then  generated  by  the  operations  of:  Reproduction  ,  Crossover  and  Mutation  , 
which  are  explained  in  greater  detail  in  Chapter  3. 

The  processes  of  reproduction,  crossover  and  mutation  are  repeated  until  an 
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entirely  new  population  has  been  generated.  The  old  generation  is  discarded,  and 
the  fitness  of  each  individual  in  the  new  generation  is  determined.  The  sum  of  the 
fitness  values  is  calculated  for  the  new  generation,  and  reproduction,  crossover  and 
mutation  are  ready  to  begin  again.  Statistics  are  maintained  to  monitor  such  things 
as  the  minimum,  maximum,  and  average  fitness  values  for  each  generation.  The 
process  is  halted  when  an  appropriate  convergence  criteria,  such  as  the  average 
fitness  value  for  a  generation,  reaxrhes  a  specified  level. 

While  the  mechanics  of  genetics  algorithms  are  fairly  easy  to  understand,  the 
underlying  reasons  for  their  success  are  more  subtle.  This  phenomenon  is  explained 
by  the  concept  of  schemata  templates  and  the  building  block  hypothesis.  Binary 
schemata  templates  are  those  in  which  some  bit  positions  are  important,  while 
other  bit  positions  are  unimportant,  or  don’t  care  bits.  For  example  the  schemata 
template  10*1**0,  may  be  thought  of  as  having  4  important  bits,  with  3  intervening 
don’t  care  bits.  For  this  schemata,  the  actual  bit  strings  1001000,  1011100,  and 
1001010  would  all  have  an  equivalent  value,  while  the  bit  string  1010110  would  have 
an  inferior  value. 

It  is  postulated  that  useful  schemata,  which  result  in  higher  fitness  values, 
concentrate  in  short  building  blocks.  These  building  blocks  are  then  sampled  and 
recombined  by  the  operations  of  reproduction  and  crossover  at  a  much  higher  fre¬ 
quency  than  those  of  inferior  combinations.  The  mutation  operation  assists  by  in¬ 
troducing  new  information  at  a  rate  that  will  not  disrupt  the  building  block  process. 
Goldberg  works  through  some  examples  to  illustrate  this  process  in  reference  [7], 
and  in  reference  [8,  page  1570]  he  asserts  that; 

Holland’s  schema  theorem  places  the  theory  of  genetic  algorithms  on 
rigorous  footing  by  calculating  a  bound  on  the  growth  of  useful  simi¬ 
larities  or  building  blocks  .  Moreover,  it  has  been  recently  conjectured. 
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that  global  convergence  of  certain  types  of  GA  is  probabilistically  poly¬ 
nomial.  If  shown  correct,  this  conjecture  could  revolutionize  the  global 
solution  of  almost  any  problem. 

1.5  Optimization  via  Genetic  Algorithm  -  An 
Example 

A  simple  optimization  problem  with  one  independent  parameter  will  now  be  used 
to  demonstrate  the  genetic  algorithm  explained  in  the  previous  section.  Assuming 
that  the  value  of  the  fitness  function  is  represented  exactly  by  the  following  equation 
in  the  parameter  x  : 

fitness  =  10.0  —  (x  —  2.5)^  +  5  sin(47rx)  (1-7) 

which  is  plotted  in  figure  1.1. 

While  simple  in  construction,  the  closely  spaced  relative  maxima  and  min¬ 
ima  in  this  example  present  a  non-trivial  optimization  problem  for  most  search 
techniques.  Using  MATLAB  to  solve  this  function  numerically  for  the  three  highest 
peaks  yields  the  following: 

X  value  peak  fitness 
2.1245  14.8589 

2.6243  14.9844 

3.1242  14.6101 

Thus  for  this  example,  the  maximum  fitness  occurs  at  the  center  peak  with  a  value 
of  14.9844,  and  the  optimum  value  of  x  is  2.6243. 

Using  a  population  size  of  20,  and  a  16  bit  string  to  represent  values  of  x 
ranging  from  0.0  to  5.0,  nine  generations  of  the  genetic  algorithm  were  performed 
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Rtness  Function 
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Figure  1.1:  Fitness  Function  for  Genetic  Algorithm  Example 

using  the  MATLAB  program  GATEST,  which  is  described  in  the  Appendix.  Fig¬ 
ure  1.2  shows  the  initial  population  of  random  values  for  generation  number  1.  The 
four  plots  in  this  figure  are  generated  for  each  generation  and  show  the  fitness  and 
parameter  values  for  each  individual  in  the  left  hand  plots,  as  well  as  maximum 
fitness  and  average  fitness  statistics  for  each  generation  in  the  right  hand  plots.  For 
the  initial  generation,  there  is  no  statistical  data  available  for  the  right  hand  plots. 

Figure  1.3  shows  the  results  after  five  generations  of  reproduction,  crossover 
and  mutation,  which  produces  an  optimum  value  for  x  of  2.6397,  correct  to  within 
.58%.  Figure  1.4  then  shows  the  results  after  nine  generations,  which  produces  an 
optimum  value  for  x  of  2.6276,  correct  to  within  .12%.  Figures  1.5  and  1.6  show 
in  greater  detail  the  history  of  the  maximum  fitness  value  and  best,  or  optimum  , 
parameter  with  increasing  generation. 

Note  that  five  generations  represent  only  100  trial  parameters,  and  nine 
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Figure  1.2:  Initial  Generation  for  Genetic  Algorithm  Example 


generations  represent  only  180  trials.  For  comparison,  a  purely  random  search 
would  require  something  on  the  order  of  5000  trials  to  achieve  an  accuracy  of  ,1% 
over  the  same  range,  and  most  gradient  search  techniques  would  have  extreme 
difficulty  optimizing  this  function. 
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Lgure  1.3: 


Results  After  Five  Generations  for  Genetic  Algorithm  Example 


Figure  1.4:  Results  After  Nine  Generations  for  Genetic  Algorithm  Example 
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Figure  1.5:  Maximum  Fitness  vs  Generation  for  Genetic  Algorithm  Example 


Genetic  Algorithm  —  Best  Parameter 


Generation 


Figure  1.6:  Best  parameter  vs  generation  for  Genetic  Algorithm  example 
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CHAPTER  2 

APPLICATION  SYSTEM  MODEL 


A  general  feedback  control  scheme  for  a  pitch  plane  missile  autopilot  controller  is 
shown  in  figure  2.1.  The  sensors  provide  feedback  variables  to  the  controller  which 
may  be  some  combination  of  state  and  response  variables.  Regarding  the  pitch 
plane  missile  dynamics,  most  actual  missile  autopilot  controllers  accept  normal 
acceleration  commands  [rje)  from  the  guidance  loop,  and  utilize  accelerometer  [rjg) 
and  rate  gyro  (0,)  feedback  information  to  determine  the  actuator  command  (^c) . 

The  actual  airframe  and  actuator  and  even  sensor  dynamics  may  be  highly 
non-linear  over  the  flight  regime  of  any  missile.  However,  over  small  ranges  of  flight 
conditions,  where  parameters  such  as  Mach  number,  dynamic  pressure,  fuel  weight, 
thrust,  etc.,  may  be  considered  constant,  linearized  models  can  be  developed  for 
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Figure  2.1:  Missile  Autopilot  Feedback  Control  Scheme 
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use  in  the  control  system  synthesis.  These  models  can  produce  controllers  that 
are  effective  within  some  range  of  the  assumed  flight  conditions.  Usually  several 
conditions  are  studied,  and  the  resulting  controller  parmeters  are  then  scheduled 
over  the  desired  flight  envelope. 

The  approach  used  in  this  thesis  will  be  to  develop  a  linearized  airframe  and 
actuator  model  and  assiime  hypothetical  values  that  correspond  to  a  single  flight 
condition.  A  state  space  model  of  this  system  will  be  developed,  with  the  states 
being: 

Xi  =  a  =  angle  of  attack 
Xi  =  6  =  pitch  angular  rate 

xs  =  6  =  actuator  deflection 

Full  state  feedback  will  be  assumed,  with  angle  of  attack  (a)  used  in  place 
of  normal  acceleration  (77)  as  the  autopilot  tracking  variable.  This  is  done  for 
simplification  in  the  development  of  the  system  model,  but  in  no  way  compromises 
the  applicability  of  the  results,  since  it  will  be  seen  that  in  the  linear  model  7/  is 
just  a  linear  combination  of  a  and  S.  For  simplicity  also,  the  sensor  dynamics  will 
be  ignored,  and  the  sensor  transfer  functions  will  assumed  to  be  unity  so  that  the 
feedback  states  are  the  actual  modelled  states. 

The  following  sections  develop  the  equations  for  the  linearized  airframe  and 
actuator  models  and  provide  the  hypothetical  flight  dependent  parameters  used  in 
the  ensuing  analysis.  A  state  space  system  model  is  then  developed,  along  with 
the  discrete  time  model  required  for  digital  simulation.  Finally,  the  measures  of 
performance  that  will  be  used  in  the  optimization  studies  are  developed  as  applied 
to  this  model. 
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2.1  Linearized  Missile  Airframe  Dynamics 

The  6  degree-of-freedom  linearized  aerodynamics  model  for  a  missile  (or  aircraft) 
rigid  body  are  developed  in  section  2.6  of  reference  [5].  We  will  extract  the  pitch 
plane  equations  from  this  development,  with  some  further  simplifications  and  changes 
in  nomenclature.  First  we  define  the  pitch  plane  dynamic  variables: 

a  -  angle  of  attack 
9  -  pitch  angle 

7  -  flight  path  angle 

V  -  velocity  vector 
T)  -  body  normal  acceleration 
6  -  actuator  deflection  angle 

X  -  body  centerline  axis 
z  -  body  normal  axis 

as  shown  in  figure  2.2. 

Beginning  with  the  three  axes  force  and  moment  equations,  reference  [5] 
proceeds  to  the  linearized  equations  for  small  perturbations  about  an  operating 
condition  by  assuming  constant  velocity  and  attitude.  Control  surfaces  and  engine 
thrust  are  also  assumed  to  be  trimmed  to  these  conditions.  Longitudinal  and  lateral 
dynamics,  which  are  normally  only  lightly  coupled,  are  assumed  uncoupled  so  that 
control  system  design  can  proceed  independently  in  each  channel.  The  resulting 
three  equations  for  the  longitudinal  (pitch)  axis  are  repeated  below,  but  with  a 
small  change  in  nomenclature  to  use  A  for  the  axial  force  (instead  of  X),  N  for  the 
normal  force  (instead  of  Z),  $  for  the  pitch  angular  rate  (instead  of  q),  and  6  for 
the  elevator  (or  actuator)  deflection  angle  (instead  of  6e)  [5,  page  45]: 


6u  =  AuAu  +  Aaa  —  gO  +  AeS 


(2.1) 
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Figure  2.2:  Missile  Pitch  Plane  Dynamic  Variables 

d  =  {NjV)Au  +  {NjV)a  +  9  +  (iVf;/V)<5  (2.2) 

0  =  MuAu  +  M^a  +  MqO  +  MeS  (2.3) 

where  the  variables  are  defined  as  follows: 

A  is  the  total  force  vector  acting  along 
the  body  x  axis  (axial  force) 

N  is  the  total  force  vector  acting  along 
the  body  z  axis  (normal  force) 

M  is  the  total  body  pitching  moment 
Au  is  the  change  in  the  velocity  vector  (V) 

6  is  the  pitch  angular  rate 
6  is  the  actuator  deflection  angle 

Note  that  in  these  equations,  partial  derivatives  of  the  force  and  moment 
terms  are  used.  In  aerodynamics,  these  forces  and  moments  are  usually  expressed 
in  terms  of  the  dimensionless  aerodynamic  coefficients,  C^,  Ca,  Cm  in  the  following 
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manner; 


N  =  QS^Cn 


(2.4) 


A  =  QS^Ca 


(2.5) 


M  =  QS^dCrr.  (2.6) 

where: 

Q  is  the  dynamic  pressure  (force/area) 

Sr  is  a  reference  area  that  is  airframe  dependent 
d  is  a  reference  length  that  is  airframe  dependent 

The  dynamic  pressiure,  Q,  is  a  function  of  air  density  and  missile  velocity 
(Q  =  lf2pV^)  ,  and  the  aerodynamic  coefficients  are  airframe  dependent,  and  may 
also  be  functions  of  parameters  such  as  angle  of  attack,  Mach  ninnber,  actuator 
deflection,  and  pitch  rate.  Thus  the  small  perturbation  derivative  approach  greatly 
simplifies  an  otherwise  highly  non-linear  problem.  Defining  the  following  partial 
derivatives  of  the  aerodynamic  coefficients  [15,  page  66]: 

C[fct  —  dC If  I  dot  (2"7) 


^Ns  —  dCiff  86 


(2.8) 


C™  =  dCJda 


(2.9) 


=  ac„/d6 


(2.10) 
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=  dCmldO  (2.11) 

Finally  the  force  and  moment  derivatives  used  in  the  linearized  longitudinal 
equations  can  be  defined  as: 


Na  =  {QSjmV){Ci,,  -  Ca) 

(2.12) 

Ng  =  {QSjmV)Cm 

(2.13) 

=  {QS,dlIy)Cma 

(2.14) 

Mg  =  {QS,dlIy)Cms 

(2.15) 

Mg  =  [QS,,d^)l[2IyV)C^g 

(2.16) 

For  this  thesis,  we  are  primarily  concerned  with  the  angular  rate  equations, 
and  will  make  the  additional  simplifying  assumptions: 

1.  velocity  is  constant  (Au  =  0) 

2.  gravity  can  be  ignored  {g  =  0) 

3.  aerodynamic  damping  is  negligible  {C^g  «  0) 

Applying  these  assumptions,  and  manipulating  V,  yields  the  following  simplified 
linear  angular  rate  equations: 


a  =  NaOc  +  0  +  NgS 


(2.17) 
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Figure  2.3:  Linearized  Airframe  Model  Block  Diagram 


e^g  =  Mad  +  M^e  +  Ms6  (2.18) 

As  stated  previously,  angle  of  attack  will  be  used  as  the  control  system 
tracking  variable  in  this  thesis.  It  should  be  noted  that  for  applications  where  the 
normal  acceleration  response  (?/)  is  required,  a  good  approximation  is  obtained  by 
assuming: 

r;  =  F'y  =  V{Nca  +  Ns6)  (2.19) 

Equations  2.17  ,  2.18  and  2.19  will  become  the  linearized  airframe  model  for  this 
thesis,  and  are  represented  in  block  diagram  form  in  figure  2.3 

Values  for  the  linearized  aerodynamic  constant  terms  are  assumed  for  pur¬ 
poses  of  analysis  in  this  thesis  that  represent  a  fictitious  but  unstable  airframe  at 
some  hypothetical  flight  condition  as  follows: 

Assumed  Flight  Condition  Constants 
A4=  64.11  iVa  =  .1803 
Mg  =  -62.34  Ng  =  .0738 
Mg  =  0  V  —  892  m/sec 
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Figure  2.4:  Linear  Actuator  Model  Block  Diagram. 

2.2  Linear  Actuator  Model 


A  simple  actuator  with  linear  response  characteristics  can  be  modelled  by  the  fol¬ 
lowing  differential  equation: 

6  =  -6)  (2.20) 

where: 

6c  is  the  commanded  actuator  angular  deflection 
6  is  the  actuator  response 
T  is  the  actuator  time  constant 

This  linear  actuator  model  is  also  represented  in  block  diagram  form  as  shown  in 
figure  2.4 

For  purposes  of  analysis  in  this  thesis,  the  following  actuator  time  constant 
will  be  used: 


T  =  .02seconds 


(2.21) 


2.3  State  Space  System  Model 

Combining  the  linearized  airframe  and  linear  actuator  models,  the  state  space  sys¬ 
tem  model  is  developed  by  making  the  following  state  variable  assignments: 

xi  =  a  —  angle  of  attack 

X2  =  0  =  pitch  angular  rate 

Xs  =  6  —  actuator  angular  deflection 
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Using  the  results  of  the  previous  two  sections,  and  making  the  appropriate 
state  variable  substitutions,  the  combined  state  equations  are: 


x’l  =  NaXi  +  X2  +  NsXz 

(2.22) 

X2  =  M„Xi  +  M5X3 

(2.23) 

Xg  =  -(l/r)x3  +  (l/r)^,. 

(2.24) 

and  the  output  equation  for  rj  is: 

y  =  UJV„Xi  +  VNsS 


(2.25) 


This  set  of  linear  differential  equations  can  be  written  in  the  familiar  matrix 
notation  as: 

X  =  Ax  +  Bu  (2.26) 


where: 


y  =  Cx  +  Du 


(2.27) 


X  is  the  3x1  state  vector 

u  is  the  1x1  control  variable  (a  scaler  in  this  case) 
A  is  the  3x3  system  matrix 
B  is  the  3x1  control  matrix 
C  is  the  1x3  state  output  connection  matrix 
D  is  the  1x1  control  output  connection  matrix 
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In  terms  of  the  previously  defined  variables,  these  matrices  for  the  state  space 
system  model  are: 


A  = 


No. 

1.0 

Ns 

Mo. 

0.0 

Ms 

0.0 

0.0 

-(1/0 

(2.28) 


B  = 


0.0 

0.0 

1/r 


(2.29) 


C  = 


VN^  0.0  VNs 


(2.30) 


D  =  [0.0] 


(2.31) 


2.4  Linear  Feedback  Controller 

The  linear  feedback  controller  provides  one  or  more  control  variables  that  are  formed 
as  a  linear  combination  of  the  system  states,  as  expressed  in  the  following  equation: 

u  =  —{kiXi  +  ktX2  +  AjsXs)  (2.32) 

or  written  in  matrix  form  as: 

u  =  -k  X  (2.33) 

The  gain  vector  k  consists  of  the  three  feedback  gains  that  are  to  be  determined 
by  the  control  system  design  or  optimization  process. 


26 


Figure  2.5:  Closed  Loop  Linear  Control  System 

This  notation  is  extended  to  tracking  systems  by  including  the  non-zero  input 
command  vector,  r,  as  follows: 

u  =  k  (r  -  x)  (2,34) 

where: 

r  is  the  3x1  input  state  command  vector 
X  is  the  3x1  state  vector 

u  is  the  1x1  control  variable  (a  scaler  in  this  case) 
k  is  the  1x3  gain  vector 

The  resulting  closed  loop  system,  when  the  linear  feedback  controller  is  applied  to 
the  linear  system  model,  is  shown  in  the  block  diagram  of  figure  2.5 

For  the  closed  loop  linear  control  system  of  figure  2.5  the  state  and  output 
equations  can  be  expressed  as: 


X  =  AcX  +  Bcr 


(2.35) 


y  =  CcX  +  Dcr 


(2.36) 


where  the  closed  loop  system  matrices  are  defined  in  terms  of  the  open  loop  system 
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matrices  and  feedback  gains  by  the  following  equations: 


Ac  =  A  -  Bk 


(2.37) 


Be  =  Bk 


(2.38) 


Cc  =  C 


(2.39) 


Dc  =  D 


(2.40) 


2.5  Discrete  Time  Simulation  Model 


Equations  2.35  through  2.40  establish  the  closed  loop  system  equations  in  contin¬ 
uous  form.  For  purposes  of  digital  simulation  and  evaluation,  however,  we  wish 
to  determine  the  corresponding  system  equations  in  discrete  form.  That  is,  given 
the  closed  loop  continuous  system  matrices,  and  a  simulation  time  interval  (dr),  we 
wish  to  determine  the  corresponding  discrete  matrices  for  the  following  equations: 


x(A:  +  1)  =  Adx(A:)  +  Bdr(A:) 


(2.41) 


y(A:)  =  Cdx(fc)  +  Ddr(fc)  (2.42) 

where  x(A:) ,  r(A:)  and  y{k)  represent  the  state,  input,  and  output  variables  at  the 
iteration,  and  x(A:  +  l)  represents  the  state  variables  at  the  (fc  +  l)‘^  iteration.  There 
are  many  ways  to  convert  from  the  continuous  to  the  discrete  system  equation,  but 
a  standard  method  recommended  in  several  of  the  references  [14,6,10],  and  used  in 
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this  thesis  involves  the  state  transition  matrix  #.  A  complete  derivation  is  given  in 
[14,  section  12.9],  with  the  following  results: 


Ad  =  $(dr) 

_  gAcdr 


=  I  +  Acdr  -h  Ac^dr72!  -h  ... 

(2.43) 

1 - 

1 

* _  O 

Be 

[idr  +  Acdr72!  +  Ac^dr^/3!  +  ...]  B* 

(2.44) 

Cd  =  Cc 

(2.45) 

Dd  =  Dc 

(2.46) 

The  equations  2.41  through  2.46  represent  the  discrete  time  simulation  model  for¬ 
mulation  for  tracking  systems  as  used  in  this  thesis.  Initial  values  for  all  states  are 
assumed  to  be  zero. 

2.6  Measures  of  Performance 

The  specific  measures  of  performance  that  will  be  used  in  forming  the  objective  of 
the  optimization  problem  in  this  thesis  are: 

1.  quadratic  performance  index 


2.  rise  time  specification 
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3.  peak  actuator  effort  specification 

4.  settling  error  specification 

Each  of  these  measures  is  generated  from  the  time  response  of  the  discrete  time 
simulation  model  to  a  unit  step  input  in  angle  of  attack.  These  measures  of  per¬ 
formance  are  used,  either  individually  or  in  weighted  combinations,  to  create  a 
single  cost  function  that  can  be  related  to  a  fitness  function  for  use  in  the  genetic 
algorithm. 

The  quadratic  performance  index  used  is  the  same  as  that  of  the  Linear 
Quadratic  Regulator  problem,  which  was  given  previously  in  continuous  form  in 
equaton  1.1.  For  the  discrete  time  simulation  model,  this  performance  index  is 
computed  from  the  system  discrete  time  response  by  the  following  equations: 

i{k)  =  [r(A:)  -  x(fc)]'  Q  [r{k)  -  x(A:)]  +  u'{k)B.u{k)  (2-47) 


where: 


J  = 


X).5dr[j(fc)  1)] 


k=0 


k  is  the  discrete  iteration  counter  (time) 
x(k)  is  the  state  variable  vector, 
r(k)  is  the  input  (tracking)  variable  vector, 
u(k)  is  the  control  variable  vector, 

Q  is  a  positive  semi-definite  weighting  matrix 
R  is  a  positive  definite  weighting  matrix 
N  is  the  number  of  discrete  iterations 
(corresponding  to  the  final  time) 


(2.48) 


30 


Values  used  for  the  state  and  controls  weighting  matrices  are: 


Q  = 


10.0  0  0 

0  0.1  0 

0  0  0.01 


(2.49) 


R  =  [0.01]  (2.50) 

For  the  unit  step  input,  the  tracking  variable  is  just  the  constant  command: 


m  = 


1.0  0  0 


(2.51) 


for  k  =  1  ...  N. 

Rise  time  specifications  vary  in  different  applications  and  are  often  defined 
differently  in  various  texts.  The  most  widely  accepted  definition  of  rise  time  is  the 
time  it  takes  for  a  signal  to  go  from  10%  to  90%  of  its  final  value  [14,  page  124]. 
This  definition  is  somewhat  cumbersome  to  compute,  however,  so  for  purposes  of 
this  thesis  the  rise  time,  indicated  by  will  represent  the  time  for  the  system 
angle  of  attack  response  (xi)  to  reach  80%  of  its  commanded  final  value. 

The  peak  actuator  effort  specification  is  simply  a  measure  of  the  peaJc  value 
of  the  actuator  response  (xj)  to  the  unit  step  input  in  angle  of  attack,  as  indicated 
by  the  following  expression: 


^max  —  max  ■{  ]x3]  }•  (2.52) 

The  settling  error  specification  is  a  measure  of  the  absolute  value  of  the 
difference  between  the  angle  of  attack  and  its  commanded  value  at  the  finite  final 
time  (r  =  Ndr).  That  is: 


e„t«e  =  |ri(iV)-Xi(JV)| 


(2.53) 
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CHAPTER  3 

GENETIC  ALGORITHM 
IMPLEMENTATION 


The  basic  form  of  the  genetic  algorithm  used  in  this  thesis  has  been  described 
previously  in  section  1.4.  This  chapter  will  now  describe  the  specific  implementation 
details  of  the  genetic  algorithm  as  applied  to  the  control  system  optimimization 
problem  of  this  thesis. 

3.1  Fitness  Function 

The  discrete  quadratic  performance  index  (J)  given  previously  in  equation  2.48 
forms  the  basis  for  constructing  a  fitness  function.  As  seen  previously,  the  genetic 
algorithm  acts  to  maximize  the  fitness  function,  while  in  fact  we  wish  to  mini¬ 
mize  the  quadratic  performance  index  (J).  We  can  accommodate  these  conflicting 
requirements  by  constructing  a  fitness  function  of  the  form: 

fitness  =  —  1  (3.1) 

where  C  is  a  cost  term  which  is  limited  by: 


0.5  <  C  <  100 


(3.2) 
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Figure  3.1:  Genetic  Algorithm  Fitness  Function  vs  Cost 

The  resulting  function  is  plotted  in  figure  3.1. 

The  factor  of  2  in  the  numerator  of  the  exponential  is  somewhat  arbitrary 
and  can  be  adjusted  to  give  the  fitness  function  a  good  gradient  in  the  region  of 
interest  for  a  specific  problem.  The  value  of  2  was  determined  after  a  few  trials  and 
was  found  to  work  well  for  the  optimization  problem  of  this  thesis.  The  constant 
value  1  is  subtracted  from  the  exponential  to  shift  the  function  to  the  origin  for 
large  cost  terms  (C  =>■  oo),  but  does  not  really  affect  the  optimization  process. 

For  a  simple  quadratic  performance  specification,  with  no  other  constraints, 
the  genetic  algorithm  fitness  function  is  formed  by  setting  the  cost  term  in  equation 
3.1  equal  to  the  quadratic  performance  index: 


C  =  J 


(3.3) 
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3.2  Performance  Penalties 


Additional  design  specification  constraints  can  be  incorporated  into  the  genetic  algo¬ 
rithm  optimization  problem  by  adding  appropriately  weighted  performance  penal¬ 
ties  to  the  cost  term  C  as  shown  in  the  following  equation: 

C  =  J  +  Pg  +  Paettle  +  Prise  {3-4) 


where: 

Pg  is  an  actuator  peak  response  penalty 
Psettie  is  a  settling  error  penalty 
Prise  is  a  rise  time  penalty 

These  performance  penalties  and  their  associated  weighting  factors  are  given 
by  the  following  expressions: 


Prise  = 


10,0  [trise  -  <f>rise)  if  ^rise  >  <f>  rise 

0.0  otherwise 


where: 


Ps  =  { 


20.0  {6,nax  -  <i>s) 

0.0 


if  >  4>S 

otherwise 


P lettle  — 


10.0  (^settle  4^settle)  if  ^settte  ^  ^settle 

0.0  otherwise 


^g  is  the  design  actuator  peak  response  specification 
<f>rise  is  the  design  rise  time  specification 
<f>settie  is  the  allowable  settling  error  specification 


(3.5) 


(3.6) 


(3.7) 
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3.3  Reproduction 

Reproduction  is  the  process  of  selecting  two  individuals  from  the  current  population 
for  mating.  Individuals  are  selected  using  a  random  function  that  gives  preferential 
weighting  to  the  individual’s  relative  fitness  within  the  population.  The  algorithm 
used  in  this  thesis  is  similar  to  the  simple  genetic  algorithm  described  in  [7]  and  is 
implemented  in  the  following  manner: 

Reproduction  Algorithm 

1.  Determine  the  sum  of  the  fitness  values  for  the  entire  population. 

2.  Rank  the  individuals  from  highest  to  lowest  fitness. 

3.  Determine  the  cumulative  fitness  values  of  each  individual,  ranked  from  high¬ 
est  to  lowest.  These  numbers  should  range  from  zero  to  the  sum  of  the  fitness 
determined  in  step  one. 

4.  Generate  a  random  number  from  a  uniform  distribution  ranging  from  zero  to 
the  sum  of  the  fitness. 

5.  Select  parent  1  for  mating  whose  cumulative  fitness  value  is  the  nearest  one 
less  than  or  equal  to  the  random  number. 

6.  Remove  parent  1  from  the  available  population  for  mating,  and  repeat  the 
steps  above  to  select  parent  2  . 

3.4  Crossover 

Crossover  is  the  process  of  combining  the  parameter  strings  from  two  parents  to 
form  new  offspring  individuals.  It  is  performed  in  a  manner  independent  of  the 
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actual  numbers  that  those  strings  represent.  For  the  multiparameter  optimization 
problem  of  this  thesis,  the  48  bit  parameter  strings  consist  of  three  16  bit  words, 
eaeh  of  which  represents  one  of  the  gain  parameters  A;i,  or  k^.  Each  parameter 
word  is  scaled  separately,  so  that  it  can  represent  a  value  between  the  kiower  and 
kupper,  to  the  maximum  resolution  allowed  by  16  bits.  That  is,  the  value  of  the  least 
significant  bit  for  each  parameter  word  is  given  by: 

kLSB  =  {kupper  “  -  l)  (3.8) 

In  the  genetic  algorithm  implemented  in  this  thesis,  two  offspring  (children) 
are  generated  from  two  parents  in  the  following  manner: 

Crossover  Algorithm 

1.  Generate  a  random  number  uniformly  distributed  between  0.0  and  1.0.  If  this 
number  is  less  than  the  predetermined  probability  of  crossover,  Pero»«>  then 
copy  the  parent  parameter  word  directly  into  the  child,  and  skip  to  the  next 
parameter  word.  Otherwise,  continue. 

2.  Generate  a  random  number  uniformly  distributed  between  1  and  the  number 
of  bits  in  a  parameter  word. 

3.  Swap  the  portion  of  the  parameter  words  between  the  two  parents  that  have 
bit  locations  higher  than  this  random  number. 

4.  Repeat  the  above  steps  for  every  parameter  word  in  the  parameter  string. 
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The  following  example  is  given  to  illustrate  this  process  for  two  single  pa¬ 
rameter  word  parent  strings  that  are  represented  by  the  following  bit  sequences: 

Pi  =  1010101010101010  (3.9) 

P2  =  1111000011110000  (3.10) 

Assuming  that  the  random  crossover  site  ntimber  is  6,  the  following  two  single 
parameter  word  child  strings  would  be  generated: 

Cl  =  101010  00  11110000  (3.11) 

C2  =  1111001010101010  (3.12) 


3.5  Mutation 

Mutation  is  a  further  randomizing  process  that  occurs  at  relatively  low  probability. 
After  reproduction  and  crossover  have  generated  an  offspring,  mutation  may  change 
the  value  of  any  bit  in  the  resulting  child  strings  from  0  to  1,  or  1  to  0.  Mutation 
can  be  implemented  in  several  ways,  but  the  straightforward  method  implemented 
in  this  thesis  is  as  follows: 

Mutation  Algorithm 

1.  For  each  bit  in  the  child  parameter  strings,  generate  a  random  number  uni¬ 
formly  distributed  beteen  0.0  and  1.0. 

2.  If  the  random  number  for  any  bit  in  the  child  parameter  string  is  less  than 
the  predetermined  probability  of  mutation,  Pm.utation'>  change  the  value  of  that 
bit.  Otherwise  leave  the  value  of  the  bit  unchanged. 
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3.6  Algorithm  Enhancements 

The  genetic  algorithm  described  in  the  previous  sections  follows  closely  the  simple 
genetic  algorithm  described  in  [7].  There  are,  however,  several  notable  differences, 
or  enhancements,  to  the  genetic  algorithm  in  this  thesis,  which  for  this  problem  at 
least,  improved  its  overall  performance.  These  enhancements  include: 

Algorithm  Improvements 

1.  Ensuring  that  no  individual  is  mated  with  itself,  if  crossover  is  to  be  performed. 

2.  Copying  the  best  2  individuals  directly  into  the  next  generation,  without  per¬ 
forming  crossover  or  mutation.  These  individuals  are  still  available  for  mating, 
but  reduce  the  total  number  of  matings  required  to  fill  the  next  generation. 
This  ensures  that  the  maximum  fitness  function  is  monotonic,  and  that  the 
best  genes  are  never  lost. 

3.  Performing  crossover  separately  on  each  parameter  word  in  the  multiparam¬ 
eter  string. 

4.  Scaling  each  parameter  word  separately  to  take  advantage  of  a  priori  infor¬ 
mation,  increase  accuracy,  and  accelerate  convergence. 

5.  Modifying  the  random  distribution  of  the  crossover  site  selection  to  favor  high 
order  bits  in  early  generations,  gradually  shifting  to  favor  low  order  bits  in 
later  generations.  The  idea  is  to  force  the  genetic  algorithm  to  make  big 
decisions  early,  then  narrow  its  focus  to  fine  tune  the  final  answer.  (This 
experimental  scheme  was  implemented  and  showed  promise,  but  was  not  used 
in  the  final  optimization  runs  presented  later  in  this  thesis.) 
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It  should  be  noted  that  several  other  algorithm  enhancements  were  attempted 
during  the  course  of  this  investigation  that  did  not  demonstrate  any  improvement 
in  the  performance  of  the  algorithm.  Some  of  these  even  degraded  its  overall  per¬ 
formance.  For  completeness,  the  algorithm  modifications  which  did  not  work  are 
listed  below: 


Unsuccessful  Algorithm  Enhancements 

1.  Modifying  the  parameter  scale  factors  as  a  function  of  generation  number, 
attempting  to  squeeze  the  upper  and  lower  limits  toward  the  optimum  values. 

2.  Saving  the  best  individuals  out  of  the  combined  old  and  new  generations  after 
each  round  of  reproduction,  crossover  and  mutation. 

3.  Generating  2  new  individuals  at  random  during  each  generation,  and  reducing 
the  number  of  individuals  generated  by  mating  by  2. 

4.  Initializing  a  portion  of  the  initial  population  to  set  strings  in  order  to  repre¬ 
sent  primitive  Walsh  function  patterns  in  the  initial  population. 


3.7  Convergence  Criteria 

In  many  optimization  problems,  it  wotild  be  desired  to  determine  a  criteria  for 
algorithm  convergence.  One  such  method,  suggested  in  reference  [8],  is  to  terminate 
the  algorithm  when 

1.  The  average  fitness  of  the  population  is  within  1%  of  the  best  fitness  for  that 
generation,  and  ... 


2.  a  minimum  of  20  generations  have  been  examined. 
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In  the  example  given  in  this  reference  for  an  LQR  controller  optimization  problem, 
it  is  observed  that  this  convergence  criteria  was  reached  after  about  65  generations. 
It  is  also  observed,  however,  that  after  about  15  generations,  the  best  fitness  is 
within  a  few  percent  of  its  final  value. 

In  the  controller  optimization  problem  of  this  thesis,  we  are  interested  in 
obtaining  the  best  fitness  value  and  are  less  concerned  with  the  average  fitness 
evolution  of  the  population  taken  as  a  whole.  For  this  reason,  as  well  as  to  conserve 
the  computer  time  required  in  order  to  produce  a  given  set  of  results,  the  genetic 
algorithm  was  executed  for  a  fixed  number  of  generations  for  each  optimization 
run.  It  was  found  that  17  generations  were  sufficient  to  produce  results  within  a 
few  percent  of  the  peak  fitness  values  for  the  cases  studied  in  this  thesis. 
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CHAPTER  4 

AUTOPILOT  CONTROLLER 
OPTIMIZATION 


This  chapter  presents  the  results  of  the  genetic  algorithm  application  to  the  missile 
autopilot  control  system  optimization  problem  posed  by  this  thesis.  First,  the  stan¬ 
dard  steady  state  LQR  solution  is  presented  and  the  resulting  closed  loop  system 
response  analyzed.  Results  for  the  genetic  algorithm  optimization  of  the  analagotis 
problem,  using  a  linear  quadratic  performance  measure,  are  then  presented  and 
compared.  Additional  design  specifications  are  next  imposed  to  meet  peak  actuator 
effort,  settling  error,  and  rise  time  constraints.  The  results  of  the  genetic  algorithm 
optimization  for  each  successive  combination  of  constraints  is  presented,  culminat¬ 
ing  in  a  controller  which  satisfies  all  of  the  constraints.  It  is  posited  that  the  resulting 
optimal  controller  is  one  for  which  analytic  solutions  are  not  available.  Finally,  the 
performance  of  the  genetic  algorithm  is  examined  by  comparing  its  solutions  to  a 
set  of  true  optimum  computer  solutions  found  with  the  help  of  the  Nelder-Meade 
numerical  search  algorithm.  The  results  of  the  genetic  algorithm  optimization  stud¬ 
ies  were  used  to  initialize  the  Nelder-Meade  algorithm,  and  a  numerical  search  was 
continued  to  find  the  true  optimum  values. 

The  computer  generated  results  presented  in  this  thesis  were  obtained  using 
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the  MATLAB  software  package,  operating  on  a  20  Mega  Hertz  386  personal  com¬ 
puter.  MATLAB  standard  functions,  as  well  as  functions  provided  in  the  Control 
Toolbox  and  Robust  Control  Toolbox,  were  utilized.  Several  MATLAB  programs 
and  numerous  special  supporting  functions  were  written  in  order  to  generate  these 
results,  and  are  provided  in  the  Appendices  for  reference.  The  MATLAB  function, 
FMINS,  was  used  to  perform  the  Nelder-Meade  numerical  optimization  studies. 


4.1  Steady  State  LQR  Solution 

The  standard  steady  state  Linear  Quadratic  Regulator  solution  to  the  continuous 
state  space  system  model,  described  in  section  2.3,  was  determined  by  solving  the 
associated  Algebraic  Riccati  equation.  This  was  performed  in  MATLAB,  using  the 
functin  LQR2  provided  with  the  MATLAB  Control  Toolbox.  The  resulting  feedback 
controller  gains  are  shown  below. 

LQR  Controller  Gains 
A:i  =  -  34.1723 
ki  =  -  3.6403 

fca  =  2.3434 

The  resulting  closed  loop  discrete  time  system,  which  incorportes  the  LQR 
feedback  gains,  was  constructed  and  analyzed  using  the  MATLAB  program  AP5. 
Analysis  of  this  analytic  controller  solution  provides  a  confidence  check  for  the 
simulation  programs  and  serves  as  a  baseline  design  for  estimating  appropriate 
ranges  of  the  gain  parameters  to  be  used  in  the  genetic  algorithm. 

The  state  variable  time  response  of  this  system  was  generated  by  inserting  a 
unit  step  input  into  the  angle  of  attack  command  variable  (ri  =  1).  The  resulting 
transient  response  tends  to  be  very  active  for  the  first  few  tenths  of  a  second. 
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then  settles  toward  the  steady  state  values  within  a  few  seconds.  For  simulation 
purposes,  a  dual  simulation  time  interval  approach  was  implemented.  First,  the 
transient  portion  of  step  response  was  generated  using  a  simulation  time  interval 
(dr)  of  .001  seconds,  then  the  state  transition  matrix  was  recomputed,  and  the 
steady  state  portion  of  the  step  response  was  generated  using  a  time  interval  of 
.01  seconds.  This  allows  accurate  integration  during  the  period  of  high  system 
dynamics,  yet  conserves  on  the  total  computer  time  required  to  evaluate  a  given 
configuration. 

Figure  4.1  shows  the  ax:tuator  time  response  during  the  transient  period  from 
0  to  .2  seconds,  from  which  the  peak  actuator  response  can  be  determined.  Figure 
4.2  shows  the  angle  of  attack  response  out  to  its  assumed  steady  state  value  at 
5.0  seconds.  From  this  angle  of  attack  history,  the  rise  time  and  settling  error 
parameters  can  be  determined. 

Bode  plots  for  the  open  loop  system,  evaluated  with  the  LQR  feedback  gains, 
are  also  obtained  in  APS  by  using  the  MATLAB  Control  Toolbox  function  BODE. 
Likewise,  gain  margin  (G„»arffin),  and  phase  margin  {<f>margin)  for  this  system  are 
computed  using  the  Toolbox  function  MARGIN.  Figure  4.3  shows  the  Bode  plots 
for  the  LQR  feedbaick  gain  open  loop  system. 

According  to  optimal  control  theory,  the  LQR  controller  design  should  ex¬ 
hibit  an  infinite  gain  margin  and  a  minimum  of  60  degrees  of  phase  margin  [1].  The 
phase  margin  can  be  obtained  from  the  Bode  plots,  as  the  phase  angle  when  the 
gain  crosses  0  dB.  Likewise,  the  gain  margin  is  obtained  from  the  Bode  plots,  as  the 
gain  value  when  the  phase  angle  crosses  180  degrees.  From  these  plots,  it  appears 
that  the  theoretical  stability  margins  have  been  achieved.  The  phase  margin  is 
somewhat  greater  than  60  degrees,  and  the  gain  margin  appears  to  be  infinite,  since 
the  phase  appears  to  approach  90  degrees  asymptotically  and  will  never  cross  180 
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Radians  /  sec 


Figure  4.3:  Bode  Plots  for  the  LQR  Controller 

degrees.  Since  no  physical  system  can  have  infinite  gain  margin,  we  will  henceforth 
indicate  only  that  this  gain  margin  exceeds  90  dB. 

Using  the  computer  to  numerically  analyze  these  time  and  frequency  re¬ 
sponse  histories,  the  following  performance  and  stability  margin  characteristics  are 
extracted  and  summarized  for  the  LQR  controller  design. 

LQR  Controller  Characteristics 


^max  — 

6.826 

.079 

trise  ~ 

.151 

^margin  ~ 

72.77  deg. 

^margin  “ 

(>  90  dB  ) 
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4.2  Genetic  Algorithm  Optimization  for  the  LQ 
Performance  Index 

The  genetic  algorithm  described  in  Chapter  3  was  implemented  using  the  MAT- 
LAB  program  APGEN3,  and  configured  to  optimize  for  the  Linear  Quadratic  per¬ 
formance  index  (J),  described  in  section  2.6.  This  was  accomplished  by  setting  the 
GA  cost  function  equal  to  the  LQ  performance  index  and  computing  the  GA  fitness 
function  as  described  in  section  3.1.  The  system  resulting  from  this  optimization 
problem  will  hence  be  referred  to  as  the:  J  Controller  .  This  is  the  genetic  algo¬ 
rithm  equivalent  of  the  LQR  design  of  the  previous  section,  using  identical  Q  and 
R  weighting  matrices.  In  this  case,  however,  numerical  analysis  of  the  discrete  time 
simulation  response  will  be  used  to  compute  and  optimize  the  GA  fitness  function. 
Differences  are  to  be  expected  due  to  the  discrete  time  simulation  intervals  (.001 
and  .01  seconds)  and  the  finite  integration  period  (5  seconds)  used  for  evaluation 
of  the  performance  index. 

Snapshots  of  the  genetic  algorithm  optimization  process  are  shown  in  fig¬ 
ures  4.4  through  4.8.  The  four  plots  in  each  of  these  figures  are  produced  after  each 
generation  of  30  individuals,  and  record  the  following  information: 

1.  The  fitness  value  for  each  individual  in  the  current  generation. 

2.  The  three  gain  parameter  values  {kijk^fks)  for  each  individual  in  the  current 


generation. 
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Generation  Number  1 


Figure  4.4:  Genetic  Algorithm  Optimization  -  Initial  Generation  -  J  Controller 


Generation  Number  5 


Figure  4.5:  Genetic  Algorithm  Optimization 


Generation  5  -  J  Controller 
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Generation  Number  9 
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Figure  4.6:  Genetic  Algorithm  Optimization  -  Generation  9  -  J  Controller 


Generation  Number  13 


Figure  4.7:  Genetic  Algorithm  Optimization  -  Generation  13  -  J  Controller 
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Generation  Number  17 


Figure  4.8:  Genetic  Algorithm  Optimization  -  Generation  17  -  J  Controller 

As  noted  previously  in  section  3.4,  each  parameter  word  in  the  genetic  algo¬ 
rithm  optimization  process  is  scaled  separately,  so  that  it  can  only  represent  gain 
values  between  kiower  ^-nd  kipper-  The  results  of  the  steady  state  LQR  solution  of 
the  previous  section  provide  a  good  indication  of  both  the  magnitude  and  sign  that 
should  be  placed  on  these  limits.  Indeed,  the  resulting  closed  loop  system  is  stable 
only  if  A:i  and  k2  are  negative,  and  ks  is  positive.  Thus  to  accelerate  the  genetic 
algorithm  convergence,  the  following  upper  and  lower  values  were  chosen  for  scaling 
the  parameter  words: 


parameter 

Slower 

k 

^uppcr 

ki 

-100.0 

0.0 

k2 

-50.0 

0.0 

ks 

0.0 

25.0 

It  should  be  noted  however,  that  it  is  not  necessary  to  limit  the  gain  param- 
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eters  to  the  correct  sign  values  in  order  for  the  genetic  algorithm  to  work.  Stability, 
as  a  hard  constraint,  is  contained  within  the  quadratic  performance  index  J.  Un¬ 
stable  systems,  due  to  an  incorrect  feedback  gain  sign,  result  in  an  extremely  large 
value  of  J,  and  are  appropriately  penalized  by  the  reproduction  function  of  the  ge¬ 
netic  algorithm.  In  trials  where  kiower  was  set  to  -100  and  kipper  to  +100  for  all 
three  gain  parameters,  the  genetic  algorithm  converged  to  the  stable  sign  region  for 
all  parameters  within  a  few  generations. 

Several  other  parameters  that  are  used  to  control  the  genetic  algorithm  pro¬ 
cess  were  determined  empirically,  and  were  set  as  follows: 


GA  Parameters 

population  size 

=  30.000 

probability  of  crossover  (pcrots) 

=  .800 

probability  of  mutation  {pmutation) 

=  .025 

Figure  4.9  shows  the  maximum  fitness  vs.  generation  in  greater  detail  for  the 
J  Controller  optimization,  and  figure  4.10  shows  best  gain  parameters  plotted  vs. 
generation.  It  is  seen  that  the  genetic  algorithm  converges  rapidly  in  this  case,  and 
is  within  a  few  percent  of  its  final  gain  values  within  3  generations.  The  feedback 
controller  gains  resulting  after  17  generations  are  shown  below, 

J  Controller  Gains 
ki  =  -  48.9113 
k2  =  -  6.2638 

ks  = 


1.6400 


Generation 


Figure  4.10:  Genetic  Algorithm  -  Best  Parameters  -  J  Controller 
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Figures  4.11  and  4.12  show  the  time  response  plots  for  the  resulting  J  Con¬ 
troller.  Analysis  of  this  controller  design  indicates  the  following  characteristics: 

J  Controller  Characteristics 


^max 

=  9.7258 

^settle 

=  .0232 

^rise 

=  .1970 

^margin 

=  48.0  deg. 

G  margin 

=  (>  90  dB  ) 

While  these  characteristics  do  not  match  exactly  those  of  the  LQR  Con¬ 
troller,  they  are  reasonably  close  considering  the  differences  in  implementing  the 
discrete  system  equation,  the  finite  integration  interval,  and  the  numerical  integra¬ 
tion  of  the  linear  quadratic  performance  index.  The  J  Controller  will  therefore  be 
used  as  a  baseline  for  subsequent  optimization  results  which  incorporate  additional 
performance  constraints. 


4.3  Peak  Actuator  Response  Constraint 

The  system  performance  characteristics  of  the  J  Controller  compare  favorably  with 
the  characteristics  of  the  LQR  Controller  design,  having  somewhat  greater  rise  time 
and  peak  actuator  response,  but  less  settling  error.  The  J  Controller  also  displays 
adequate  robustness  characteristics,  with  gain  and  phase  margins  approaching  the 
theoretical  values  of  the  LQR  Controller.  Let  us  now  assume  that  the  actuator 
system  we  wish  to  use  has  a  physical  limit  of  ±  5  degrees  of  travel,  and  impose  this 
as  a  constraint  on  the  peak  actuator  reponse.  This  constraint  is  incorporated  into 
the  genetic  algorithm  optimization  process  by  setting  the  actuator  peak  response 
specification  term: 


<f>s  =  5.0 


(4.1) 
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Figure  4.13:  Genetic  Algorithm  -  Maximum  Fitness  -  JPg  Controller 

and  computing  the  actuator  peak  response  penalty  (/’;)  as  defined  in  equation  3.6. 
This  penalty  is  then  added  to  the  GA  cost  function,  and  the  resulting  genetically 
optimized  system  becomes  the  JPg  Controller. 

Figure  4.13  shows  the  maximum  fitness  vs.  generation  in  greater  detail 
for  the  JPg  Controller  optimization,  and  figure  4.14  shows  best  gain  parameters 
plotted  vs.  generation.  The  feedback  controller  gains  resulting  after  17  generations 
are  shown  below. 

JPg  Controller  Gains 
fci  =  -  66.0487 
ki  =  -  14.9935 

kz  — 


9.1844 
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Beat  Parameters  vs.  Generation 


Generation 


Figure  4.14:  Genetic  Algorithm  -  Best  Parameters  -  JPg  Controller 


Figures  4.15  and  4.16  show  the  time  response  plots  for  the  JPg  Controller. 
Analysis  of  this  controller  design  indicates  the  following  characteristics: 

Jpg  Controller  Characteristics 
^maz  =  4.9977 

^atUU  =  .1116 

triat  =  .3200 

^margin  83.72  deg. 

^margin  ~  90  dB  ) 

It  is  observed  that  the  peak  actuator  response  constraint  has  been  met  by 
the  new  controller  design,  and  the  stability  margin  data  show  that  robustness  has 
not  been  compromised.  However,  it  is  noted  that  the  settling  error  has  increzised 
dramatically  with  this  design.  Rising  from  approximately  2%  for  the  previous  con¬ 
troller,  the  settling  error  for  this  design  is  now  over  11%  of  the  commanded  value. 
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4.4  Settling  Error  Constraint 

Let  us  now  assume  that  the  settling  error  observed  in  the  J  Pg  Controller  is  unsatis¬ 
factory  in  terms  of  system  performance,  and  impose  this  as  an  additional  constraint 
on  the  optimization  problem.  This  constraint  is  incorporated  into  the  genetic  algo¬ 
rithm  by  adding  the  weighted  settling  error  penalty  {Psettie)  as  defined  in  equation 
3.7  to  the  GA  cost  function.  The  resulting  genetically  optimized  system,  which 
incorporates  peaJi  actuator  response  and  settling  error  constraints,  will  be  referred 
to  as  the  JPgPaettie  Controller. 

Figure  4.17  shows  the  maximum  fitness  vs.  generation  in  greater  detail  for 
the  JPsPaettie  Controller  optimization,  and  figure  4.18  shows  best  gain  parameters 
plotted  vs.  generation.  The  feedback  controller  gains  resulting  after  17  generations 
are  shown  below. 

JPsPiettie  Controller  Gains 
ki  =  -  24.7425 
fcj  =  -  7.0481 

kz  =  1.6289 

Figures  4.19  and  4.20  show  the  time  response  plots  for  the  JPsPaettu  Con¬ 
troller.  Analysis  of  this  controller  design  indicates  the  following  characteristics: 


^6^ Bettle 

Controller  Characteristics 

^max 

=  4.7937 

^settle 

=  .0377 

^rise 

=  .4400 

^margin 

=  46.7  deg. 

GffiQrgin 

=  (>  90  dB  ) 

It  is  observed  that  the  settling  error  heis  been  reduced  to  less  than  4%  by 
the  new  controller  design,  and  peak  actuator  response  constraint  has  been  simul- 


Best  Parameters:  k1»  k2,  k3 


Generation 


Figure  4.17:  Genetic  Algorithm  -  Maximum  Fitness  -  JPsP$ettu  Controller 


Best  Parameters  vs.  Generation 


Generation 


Figure  4.18:  Genetic  Algorithm  -  Best  Parameters  -  JPsPaettu  Controller 
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Figure  4.19:  Actuator  Response  for  the  JPsPtettie  Controller 


taneously  satisfied.  Stability  margin  data  also  indicate  that  robustness  has  not 
been  compromised.  However,  it  is  noted  that  the  rise  time  of  the  step  response  has 
increased  dramatically  from  .32  seconds  for  the  previous  design  to  a  sluggish  .44 
seconds  for  this  design. 

4.5  Rise  Time  Constraint 

Let  us  now  assume  that  the  sluggish  rise  time  observed  in  the  JPsPeettu  Controller  is 
unsatisfactory,  and  that  the  system  performance  specification  requires  that  the  rise 
time  be  no  greater  that  .3  seconds.  This  final  additional  constraint  is  incorporated 
into  the  genetic  algorithm  optimization  process,  by  setting  the  rise  time  specification 
term: 


4^rise  —  -3 


(4.2) 


59 


Figure  4.20:  Angle  of  Attack  Response  for  the  JPsPsettu  Controller 

and  computing  the  rise  time  penalty  {PrUt)  as  defined  in  equation  3.5.  This  penalty 
is  then  added  to  the  GA  cost  function,  and  the  resulting  genetically  optimized 
system  becomes  the  J PsPaettu Prise  Controller. 

Figure  4.21  shows  the  maximum  fitness  vs.  generation  in  greater  detail 
for  the  JPsPsettie Prise  Controller  optimization,  and  figure  4.22  shows  best  gain  pa¬ 
rameters  plotted  vs.  generation.  The  feedback  controller  gains  resulting  after  17 
generations  are  shown  below. 

Pf  Psettie  Prise  Controller  Gains 
ki  =  -  22.3987 
ki  =  -  4.5464 

kz  = 


1.6514 


Best  Parameters:  k1,  k2,  k3 


Generation 


Figure  4.21:  Genetic  Algorithm  -  Maximum  Fitness  -  J Ps P»ettiePri»e  Controller 


Best  Parameters  vs.  Generation 


Figure  4.22:  Genetic  Algorithm  -  Best  Parameters  -  J PsPaettuPrUe  Controller 
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Figures  4.23  and  4.24  show  the  time  response  plots  for  the  J PsPsettuPriee 


Controller.  Analysis  of  this  controller  design  indicates  the  following  charzicteristics: 
JPsPsettie  Prise  Controller  Characteristics 


^settle 


4.8502 

.0750 

.3000 


^margin  —  56.49  deg. 

G  margin  ~  (>  90  dB  ) 

It  is  observed  that  the  rise  time  has  been  reduced  to  .3  seconds,  which  meets 
exactly  the  system  performance  specification.  Also,  the  peak  actuator  response  is 
within  its  physical  limit  of  ±  5  degrees  of  travel.  Meanwhile,  the  settling  error 
has  increased  to  7.5%,  but  is  still  lower  than  the  11%  obtained  with  no  settling 
error  penalty,  and  may  be  acceptable.  Adequate  gain  and  phase  margins  are  also 
maintained  for  robustness  of  this  design. 


4.6  Genetic  Algorithm  Performance 

In  order  to  determine  the  accuracy  and  performance  of  the  genetic  algorithm  in 
this  application,  the  MATLAB  program  APNMOPT  was  constructed  to  find  the 
true  optimum  solution  in  each  of  the  controller  specifications  presented.  The  APN¬ 
MOPT  program  uses  the  standard  MATLAB  function,  FMINS,  which  implements 
the  Nelder-Meade  simplex  algorithm  to  find  the  minimum  of  a  multivariable  func¬ 
tion.  Since  it  is  known  that  this  algorithm  is  not  well  suited  for  finding  global 
optimums,  the  optimum  gain  values  found  previously  by  the  genetic  algorithm  were 
used  cLS  the  initial  guess  required  by  the  Nelder-Meade  algorithm. 

The  true  optimum  gain  and  fitness  function  values  obtained  by  the  Nelder- 
Meade  algorithm  for  each  of  the  four  controller  specifications  studie,  are  summarized 
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in  the  table  below.  The  genetic  algorithm  results  are  also  summarized  in  a  similar 
table  for  comparison. 

Nelder-Meade  Optimization  Summary 


Controller 

ki 

ki 

kz 

Fitness 

J 

-27.0995 

-3.5254 

0.4122 

3.2230 

JPs 

-65.1220 

-16.0587 

8.8373 

1.5807 

J  PsP settle 

-21.4221 

-6.9169 

0.7232 

1.7239 

J PsP settle  Prise 

-22.2218 

-4.4455 

1.4819 

1.3316 

GA 

Optimization  Summary 

Controller 

ki 

ki 

Fitness 

J 

-48.9113 

-6.2638 

1.6400 

3.1183 

JPs 

-66.0487 

-14.9935 

9.1844 

1.5429 

J  PsPsettle 

-24.7425 

-7.0481 

1.6289 

1.4465 

J  PsP settle  P rise 

-22.3987 

-4.5464 

1.6514 

1.2367 

The  resulting  errors  in  the  genetic  algorithm  optimized  fitness  functions  are  com¬ 
pared  with  the  true  optimum  values  obtained  with  the  Nelder-Meade  algorithm  in 
the  table  below; 


Peak  Fitness  Summary 


Controller 

GA  Optimum 

True  Optimum 

Error 

J 

3.1183 

3.2230 

3.2% 

JPs 

1.5429 

1.5807 

2.4% 

J  PsP settle 

1.4465 

1.7239 

16.1% 

J  PsP settle  P rise 

1.2367 

1.3316 

7.1% 

From  the  tables  above,  it  is  observed  that  the  genetic  algorithm  optimum 
values  are  close  to  the  true  optimum  values,  with  a  few  possible  exceptions.  For 
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the  genetically  optimized  J  Controller,  the  fej  gain  differed  considerably  from  the 
true  optimum  value.  Nevertheless,  this  seemed  to  have  little  impact  on  the  fitness 
function,  which  was  within  3.2%  of  the  true  optimum  value.  For  the  JPsPiettu  Con¬ 
troller,  the  GA  fitness  function  was  in  error  by  16.1%,  yet  the  gains  ki,  k2,  and  k^ 
were  reasonably  close  to  their  true  optimum  values. 

The  performance  characteristics  of  the  four  genetically  optimized  controllers 
have  been  previously  analyzed,  and  are  summarized  for  reference  in  the  table  below. 
Likewise,  the  performance  characteristics  of  the  controllers  using  the  Nelder-Meade 
true  optimum  gain  values  are  also  tabulated  below.  From  this  comparison,  it  is 
apparent  that  the  GA  optimized  controllers,  although  somewhat  suboptimal,  have 
good  performance  characteristics  that  refiect  the  embedded  combinations  of  design 
criteria. 


GA  Optimized  Controller  Characteristics 


Controller 

^max 

^settle 

^rise 

^margin 

J 

9.7258 

.0232 

.1970 

48.0  deg 

JPs 

4.9977 

.1116 

.3200 

83.7  deg 

J  ^ S  ^ settle 

4.7937 

.0377 

.4400 

46.7  deg 

J  PsP settle  P rise 

4.8502 

.0750 

.3000 

56.5  deg 

True  Optimum  Controller  Characteristics 


Controller 

^max 

^settle 

^rise 

^margin 

J 

8.2360 

.0207 

.2000 

32.9  deg 

JPs 

5.0000 

.1015 

.3500 

82.1  deg 

J  Pd  P settle 

4.9156 

.0000 

.5200 

31.6  deg 

J  PsP settle  P rise 

5.0000 

.0679 

.3000 

53.73  deg 
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For  better  understanding  of  the  sensitivity  of  the  fitness  function  to  the 
independent  parameters,  three  comparison  plots  have  been  made  for  each  controller 
specification.  Since  the  resulting  four  dimensional  performance  surfaces  are  difficult 
to  visualize  and  impossible  to  represent  in  two  dimensions,  the  following  procedure 
was  used: 

1.  Evaluate  and  plot  the  fitness  function  vs.  ki,  while  holding  constant  the  GA 
optimum  values  for  Aij  and  k^. 

2.  Evaluate  and  overlay  on  this  plot  the  fitness  function  vs  ki,  for  constant  true 
optimum  values  of  k2  and  k^. 

3.  Repeat  1.  and  2.  above,  but  generate  plots  vs.  ^2,  while  using  constant  values 
of  ki  and  k^. 

4.  Repeat  1.  and  2.  above,  but  generate  plots  vs.  kz,  while  using  constant  values 
of  ki  and  k^. 

The  resulting  parametric  fitness  function  comparisons  are  shown  in  figures 
4.25  through  4.27  for  the  J  Controller,  figures  4.28  through  4.30  for  the  JPs  Con¬ 
troller,  figures  4.31  through  4.33  for  the  JPsPsettu  Controller,  and  figures  4.34 
through  4.36  for  the  J PsPaettu Prise  Controller.  These  figures  were  generated  us¬ 
ing  MATLAB  program  APOPLOT  after  the  GA  and  true  optimum  parameters 
were  known.  The  (*)  on  each  plot  marks  the  peak  value  for  the  true  optimum  of 
the  fitness  function  found  by  the  Nelder-Meade  algorithm.  The  (X)  on  each  plot 
marks  the  optimum  fitness  value  found  by  the  genetic  algorithm  for  that  case. 

As  indicated  by  these  results,  the  genetic  algorithm  converged  to  near  opti¬ 
mal  controller  designs,  within  17  generations,  for  each  progressively  more  difficult 
combination  of  design  constraints.  It  should  be  noted  that  the  17  generations  of 
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Figure  4.25:  Fitness  Function  vs  ki  for  the  J  Controller  Specification 

30  individuals  represent  only  478  combinations  (or  less)  of  the  three  gain  parame¬ 
ters  (allowing  for  the  fact  that  the  two  best  individuals  are  duplicated  exactly  in 
each  following  generation).  For  practical  applications,  however,  more  than  17  gen¬ 
erations  may  be  required  to  assure  accurate  convergence  of  the  genetic  algorithm. 
The  previously  cited  example  by  Goldberg  required  approximately  65  generations 
to  converge  [8].  Perhaps  a  useful  strategy  would  be  to  use  the  genetic  algorithm  for 
some  number  of  generations  to  assure  global  optimality,  then  switch  to  a  more  or¬ 
thodox  numerical  search  method,  such  as  the  Nelder-Meade  algorithm,  to  converge 
on  the  final  optimum. 


k2  gain 


Figure  4.26:  Fitness  Function  vs  k2  for  the  J  Controller  Specification 


Fitness  Function  vs  k3 


k3  gain 


Figure  4.27:  Fitness  Function  vs  fca  for  the  J  Controller  Specification 


k2  gain 


Figure  4.29:  Fitness  Function  vs  k2  for  the  JPg  Controller  Specification 


Rtness  hri  Rtness 
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Fitness  ^  Rtness 


1.8 


ire  4.34:  Fitness  Function  vs  ki  for  the  JPsPtettuPrite  Controller  Specification 


Rtness  Function  vs  k2 


Figure  4.35:  Fitness  Function  vs  kt  for  the  J Ps  P,ettie  Prise  Controller  Specification 
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CHAPTER  5 
CONCLUSIONS 


In  this  thesis  a  control  system  optimization  problem  was  formulated  with  respect  to 
a  simple  pitch  plane  missile  autopilot.  Linearized  missile  airframe  and  actuator  dy¬ 
namics  models  were  developed,  and  a  linear  state  feedback  controller  was  assumed. 
A  set  of  hypothetical  flight  parameters  were  assumed,  and  a  discrete  time  state 
space  simulation  model  for  the  closed  loop  system  was  developed. 

A  discretized  linear  quadratic  performance  measure  was  developed,  which 
can  be  computed  numerically  from  the  closed  loop  simulation  model.  Additional 
performance  specifications  were  then  imposed  to  provide  constraints  on  peak  actu¬ 
ator  response,  settling  error,  and  rise  time  performance.  It  was  posited  that  the 
resulting  multicriterion  optimization  problem  was  one  for  which  analytic  solutions 
are  not  available. 

The  genetic  algorithm  was  then  introduced,  and  a  suitable  fitness  measure 
was  developed  that  could  incorporate  the  stated  performance  measure  and  auxiliary 
constraints.  Computer  simulation  results  were  then  presented  that  demonstrate  the 
genetic  algorithm  provides  good  convergence  to  near  optimal  controller  designs  for 
successive  combinations  of  constraints.  The  final  genetically  optimized  controller 
was  designed  to  minimize  a  linear  quadratic  performance  index,  while  simultane¬ 
ously  minimizing  settling  error  and  satisfying  hard  constraints  on  peak  actuator 
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response  and  mininiuni  rise  time  performance.  The  demonstrated  success  of  the 
genetic  algorithm  in  this  application  satisfies  the  objective  of  this  thesis. 

5.1  Future  Applications 

From  this  vantage  point,  future  applications  for  the  genetic  algorithm  in  control 
system  optimization  seem  promising  and  virtually  unlimited.  A  simple  application 
system  model  was  used  in  this  thesis  primarily  to  meet  time  and  computer  resource 
constraints.  Given  sufficient  processing  resources,  however,  the  closed  loop  system 
model  evaluated  by  the  genetic  algorithm  for  this  application  could  be  replaced  by 
a  full  scale  simulation  of  arbitrary  complexity. 

Further  investigations  into  suitable  performance  measures  and  the  incorpora¬ 
tion  of  many  additional  design  contraints  are  also  warranted.  Robustness  contraints 
were  not  addressed  directly  in  this  thesisj  rather  robust  stability  margins  were 
achieved  as  a  property  of  full  state  feedback  and  application  of  the  linear  quadratic 
performance  index.  A  viable  area  of  study  would  be  to  incorporate  multivariable 
robust  control  specifications  into  the  genetic  algorithm  fitness  function. 

The  linear  feedback  controller  form  used  in  this  thesis  was  also  selected  as  a 
comfortable  starting  place.  Most  real  world  control  system  design  problems  do  not 
enjoy  the  luxury  of  full  state  feedback,  free  of  sensor  noise  and  time  delay.  Applica¬ 
tion  of  the  genetic  algorithm  optimization  approach  could  easily  be  applied  to  other 
forms  of  controllers.  Realistic  control  systems  applications  involving  reduced  order 
observers  and  optimal  stochastic  observers  should  be  investigated.  The  GA  could 
potentially  be  used  to  optimize  both  controller  and  observer  in  such  configurations. 

Finally,  there  is  no  reason  to  limit  the  genetic  algorithm  optimization  process 
to  linear  systems,  or  even  to  time  invariant  systems.  For  time  varying  systems,  it 
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may  one  day  be  possible  to  implement  a  genetic  algorithm  optimization  loop  in  real 
time.  This  depends  of  course  upon  such  things  as  the  processing  power  available, 
the  complexity  of  the  problem,  and  the  real  time  available  for  the  optimization 
loop.  It  may  be  seen  however,  that  the  genetic  algorithm  has  an  inherent  advan¬ 
tage  for  real  time  implementation  within  a  parallel  processing  environment.  This 
natural  parallelism  of  the  genetic  algorithm  can  be  exploited,  for  example,  by  in¬ 
voking  simultaneous  evaluation  of  the  fitness  function  for  all  the  individuals  (32, 
or  64,  or  128,  ...)  in  a  given  generation,  each  assigned  to  a  separate  processor. 
These  evaluations  are  the  most  time  consuming  operations  of  the  genetic  algorithm 
optimization  process.  The  serial  operations  required  to  select  and  mate  individuals 
between  generations  represent  only  a  small  fraction  of  the  total  processing  time 
requirements! 

Not  every  controller  optimization  problem  is  feasible,  but  for  those  that  are 
feasible,  the  genetic  algorithm  can  provide  a  useful  tool  for  automating  the  controller 
design  process. 
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APPENDIX  A 


MATLAB  Programs 


Several  programs  were  written  in  MATLAB  to  produce  the  analysis  and  generate 
the  various  plots  presented  in  this  thesis.  This  appendix  briefly  describes  the  main 
programs  that  were  used  and  provides  a  listing  of  these  progranos.  The  supporting 
functions  that  were  created  in  MATLAB  especially  for  these  programs  are  described 
in  the  following  appendix.  In  addition  these  programs  require  numerous  standard 
MATLAB  functions,  as  well  as  supporting  functions  from  the  MATLAB  Control 
System  Toolbox  [13,10]. 

A.l  APGEN3 

This  MATLAB  program  implements  a  genetic  algorithm  in  order  to  optimize  a  mis¬ 
sile  linear  feedback  controller.  It  provides  continuous  on  screen  plots  demonstrating 
algorithm  progress,  and  provides  automatic  or  manually  controlled  options  to  send 
output  to  a  printer.  Several  experimental  options  are  controlled  by  use  of  option 
flags,  as  explained  in  the  comments. 

A  listing  of  this  program  follows. 
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echo  on 

APGEN3  -  Autopilot  Genetic  Algorithm  Design  Program  2 

-  optimize  three  parameters  kl  fc  k2  &  k3 

-  airframe  model  only 

-  with  linear  feedback  controller 

-  numerical  cost  function  computation 
for  optimization  studies 

-  using  LQR  cost  function 

-  EVALP0P3  includes  rise  time  penalty  (  1/26/93  ) 

-  corrected  N_delta  (  2/8/93  ) 

-  Actuator  model  (  2/8/93  ) 

-  EVALP0P4  to  include  hang  off  error  penalty  (2/8/93) 

-  EVALP0P6  to  use  3  parameterr  (  2/9/93  ) 

-  Limit  squeeze  algorithm  (  2/10/93  ) 

-  Option  to  vary  crossover  site  distribution  with 
generation  number  (  2/16/93  ) 

-  EVALP0P6  to  include  actuator  max  response  penalty 
(  2/23/93  ) 

-  Add  hoe_thresh  to  EVALP0P6  (  2/24/93  ) 

-  EVALP0P7  and  STEP2S  to  use  split  dt’s  for  evaluating 
step  function  -  to  correct  actuator  norm  (  3/1/93  ) 

by  R.  Hull.  2/9/93 

This  program  implements  a  genetic  algorithm  to  optimize  the 

autopilot  controller  for  a  simplified  missile  autopilot  system. 


Special  Constants  ; 

HztR  =  2*pi:  %  Converts  Hz  to  Radians  per  Second 

walsh.string  =  [1010101010101010 
1100110011001100 
1111000011110000 
111  1111100000000 
1111111111111111 
0101010101010101 
0011001100110011 
00001  1  1100001111 
0000000011111111 
000000000000000  0] 

%  Global  Variables  : 
global  cont_plot  ; 

cont_plot  =  1  :  %  Continuous  plot  flag 

auto_plot  =  1  :  %  Automatic  plot  flag  (  to  printer  ) 

combine_flag  =  0  ;  %  Option  to  combine  old  and  new  populations 

%  and  save  the  best  individuals  of  both. 
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xover_flag  =  0  ;  %  Use  crossover  site  distribution  as  a 

%  function  of  generation  number. 

squeeze.lim  =  0  ;  %  Squeeze  limits  flag 

1,  Genetic  Algorithm  Structures 

maxgen  =  17  ;  %  total  number  of  generations 

print_gen  =  4  ;  %  number  of  generations  per  output  loop 

popsize  =  30  :  %  number  of  individules  in  population 

%  -  must  be  even  ! 

topsize  =  2  ;  %  number  (even)  of  top  ranked  individuals  to 

%  propagate  exactly  in  the  next  generation 

if  combine_flag  ==  1  %  Make  sure  topsize  is  0  if  combine  option 

topsize  =  0  :  %  is  used, 

end  ; 

newsize  =  0  ;  %  number  (  even  )  of  new  random  individuals 

%  to  include  in  each  new  generation 

numparams  =  3  ;  %  number  of  parameters  to  optimize 

paramlength  =  16  ;  %  number  of  bits  per  parameter 

stringlength  =  numparams  *  paramlength  ;  %  genetic  string  length 

%  (  bits  ) 

dt  =  .01  ;  %  Computaion  interval  -  sec. 

t.final  =  5.0  :  %  cost  &  fitness  function  final  time  -  sec. 

dts  =  .001  :  %  Small  computaion  interval  -  sec. 

ts.final  =  0.2  :  %  Small  dt  response  final  time  -  sec. 

dtl  =  .01  :  %  Small  computaion  interval  -  sec. 

tl_final  =  5.0  :  1,  Small  dt  response  final  time  -  sec. 


%  J  Controller  Specification  : 

hoe_thresh  =  100.0  ;  %  settling  error  penalty  threshold 
phi_rise  =  100.0  ;  %  rise  time  penalty  threshold 

act_thresh  =  100.0  ;  %  maximum  actuator  response  specification 

%  penalty  threshold 

%  J  Pd  Controller  Specification  : 

hoe_thresh  =  100.0  ;  %  settling  error  penalty  threshold 
phi_rise  =  100.0  ;  %  rise  time  penalty  threshold 

act_thresh  =  5.0  ;  %  maximum  actuator  response  specification 

%  penalty  threshold 
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%  J  Pd  Psettle  Controller  Specification  : 

hoe_thresh  =  0.0  ;  %  settling  error  penalty  threshold 

phi_rise  =  100.0  ;  %  rise  time  penalty  threshold 

act~thresh  =  5.0  ;  %  maiximum  actuator  response  specification 

%  penalty  threshold 


%  J  Pd  Psettle  Prise  Controller  Specification  : 

hoe_thresh  =0.0  %  settling  error  penalty  threshold 

phl_rise  =  .30  %  rise  time  penalty  threshold 

act_thresh  =6.0  %  maximum  actuator  response  specification 

y,  penalty  threshold 

%  Verify  controller  specification 
pause 

pcross  =  .8  ;  %  probability  of  crossover 

pmutation  =  .025  %  probability  of  mutation 

kg.scale  =  [  -100  0 

-60  0 

0  25  ]  ;  X  gain  parameter  scales 

init_scale  =  kg.scale  %  used  for  limit  squeeze  algorithm 
squeeze_factor  =  0.9  %  used  for  limit  squeeze  algorithm  (<  1.0) 

squeeze.buffer  =2.0  %  used  for  limit  squeeze  algorithm  (  >  0  ) 


%  Flight  Condition  Constants  : 

M_alpha  =64.11; 

M_delta  =  -62.34; 

M_thetadot  =  0; 

N.alpha  =  .1803; 

N_delta  =  .0738; 

Tau.act  =  .02  ;  %  Actuator  first  order  time  constant 

Velocity  =  892.0;  %  meters/sec 

%  Units  Conversions 


y.  Define  Airframe  +  First  Order  Actuator  State  Space  Realization  : 

A  =  [  -  N.alpha  1.0  -  N_delta 

M_alpha  0  M_delta 

0  0  -  (  1  /  Tau_act  )  ]  ; 

B  =  [  0 

0 

(  1  /  Tau_act  )  ] 


C  =  [  (  N^alpha  *  Velocity  )  0  (  N_delta  *  Velocity  )  ]  ; 


D  =  [  0  ]  ; 

%  Determine  LQR  design  for  basic  plant  with  first  order  actuator 

qmatrix  =  [  10  0  0 

0.10 
0  0  .01  ]; 

rmatrix  .01; 

%  Initialize  time  values 
t  =  0:dt:t_final  ; 
ts  =  0:dts:ts_final  ; 
tl  =  ts_f inal:dtl:tl_f inal  ; 

%  Set  screen  display  format  to  suppress  excessive  line  feeds 
format  compact  ; 

%  Experimental  Option  to  Initialize  Population  using  Walsh  Strings 
%  new.chrom  =  [  ]  ; 

%  for  iw  =  1  :  popsize 
%  new_chrom  =  [  new^chrom  ;  . . 

%  walsh_string(  iw,  :  )  walsh_string(  iw,  ;  )  walsh_string(  iw, 
y.  pause 
%  end  ; 

y«  Initialize  New  Population  to  Random  Strings 

new.chrom  =  fix(  1.9999999  ♦  rand(  popsize,  stringlength  )  ) 

%  Clear  generational  statistics 

gen_avg_f itness  =  [  ]  ; 

gen_max_f itness  =  [  ]  ; 

gen_max_param  =  [  ]  ; 
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geii_max_chrom  =  [  ]  ; 

gen_inin_litness  =  [  ]  : 
gen_miii_para]n  =  [  ]  ; 
geii_min_chrom  =  [  ]  : 
gen_kg_lower  =  [  ]  : 
gen_kg_upper  =  [  ]  ; 
gen  =  1  :  %  generation  number 

%  Evaluate  Fitness  of  Initial  Population 

[  new_fitness,  new_param  ]  =  evalpop7(  A,  B,  C,  D,  ts,  dts,  tl,  dtl, 
iioe.thresh,  phi.rise,  act.thresh,  popsize,  kg_scale,  new.chrom  ) 

y,  Evaluate  Initial  Population  Statistics 
avg_fitness  =  sum(  new_fitness  )  /  popsize  ; 

[  max_fitness,  imax  ]  =  max(  new_fitness  )  ; 
max_param  =  new_param(  imax,  :  )  ; 

max_chrom  =  new_chrom(  imax,  :  )  ; 

[  min.fitness,  imin  ]  =  min(  new.fitness  )  ; 
min.param  =  new_param(  imin,  :  )  ; 

min.chrom  =  new_chrom(  imin,  ;  )  ; 

%  Save  initial  generation  data  ; 

gen_avg_f itness  =  [  gen_avg_f itness  ;  avg.fitness  ]  ; 
gen_max_fitness  «=  [  gen_max_f itness  ;  max.fitness  ]  ; 
gen_max_param  =  [  gen_max_param  ;  max.param  ]  ; 

gen_max_chrom  =  [  gen_max_chrom  ;  max.chrom  ]  ; 

gen_min_f itness  =  [  gen_min_f itness  ;  min.fitness  ]  ; 
gen.min.param  =  [  gen.min.param  ;  min.param  ]  ; 

gen.min.chrom  =  [  gen.min.chrom  ;  min.chrom  ]  ; 

gen.kg.lower  =  [  gen.kg.lower  ;  kg_scale(  ;,  1  ).’  ]  : 
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geii_kg_upper  =  [  geii_kg_upper  ;  kg_scale(  2  ).’  ]  ; 

%  Display  or  plot  values  to  screen 
if  cont_plot  ==  1 

clg  : 

%  Plot  of  fitness  vs  individual  for  this  generation  : 

subplot (221)  ; 

ind_num  =  1  :  popsize  ; 

plot(  ind_num,  new_fitness  ), 

title( [’Generation  Number  int2str(gen)]  ),. 
xlabel ( ’ Individual ’ ) , . . 
ylabeK 'Fitness') , . . 
grid 


%  Plot  parameter  kl  vs  individual  for  this  generation  : 

subplot (223)  ; 

ind_num  =  1  :  popsize  ; 

axis([  0  popsize  kg_scale(l , 1)  kg_scale(l ,2)  ]  )  ; 
plot(  ind_num,  new_param( : , 1)  ), 

title ([’Generation  Number  ’,  int2str(gen)]  ),. 
xlabel ( ’ Individual ’ ) , . . 
ylabel ( ’ Parameter ;  kl ’ ) , . . 
grid 

%  Plot  parameter  k2  vs  individual  for  this  generation  : 

subplot (222)  ; 

ind_num  =  1  :  popsize  ; 

axis([  0  popsize  kg_scale(2,l)  kg_scale(2,2)  ]  )  ; 
plot(  ind_num,  new_param( ; ,2)  ), 

title( [’Generation  Number  ’,  int2str(gen)]  ),. 
xlabel ( ’ Individual ’ ) , . . 
ylabel ( ’ Parameter :  k2 ’ ) , . . 
grid 

%  Plot  parameter  k3  vs  individual  for  this  generation  : 

subplot (224)  ; 

ind_num  =  1  :  popsize  ; 

axis([  0  popsize  kg_scale(3,l)  kg_scale(3,2)  ]  )  : 
plot(  ind_num,  new_param( : ,3)  )  , 

title ( [’Generation  Number  ’,  int2str(gen)]  ),. 
xlabel ( ’ Individual ’ ) , . . 
ylabel ( ’ Parameter :  k3 ’ ) , . . 
grid 


if  auto_plot  ==  1 ,  print ,  shg ,  end  ; 
else 


generation  =  [  gen,  avg_fitness,  max_fitness,  min_fitness  ] 
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end  : 


while  1  %  operator  control  loop 

%  Genetic  Algorithm  Optimization  Loop 
lor  igen  =  1  :  print „gen 
gen  =  gen  +  1  ; 

%  Save  Old  Generation 
old_fitness  =  new^fitness  ; 
old_chrom  =  new^chrom  ; 
old_param  =  new_param  ; 
new^chrom  =  [  ]  ; 

snm.fitness  =  snm(  old^litness  )  ; 

%  Squeeze  Limits  Algorithm 

if  squeeze_lim  ==  1 

for  np  =  1  :  numparams 

last.lower  =  kg_scale(  np,  1  )  ; 

last^upper  =  kg^scale(  np,  2  )  ; 

data^lower  =  min(  old_param(  np  ))  ; 

data„upper  =  max(  old„param(  np  ))  ; 

[  best_f itness,  imax  ]  =  max(  old^fitness  )  ; 
best^param  =  old_param(  imeix,  np  )  ; 

%  Squeeze  limits  toward  best  parameter  value 

new.lower  =  best.param  -  . , 

squeeze^! actor  *  abs(  last_lower  -  best_param  ) 

new^upper  =  best_param  +  . . 

squeeze^factor  ♦  abs(  last^upper  -  best_param  ) 

%  New  limits  cannot  exclude  old  data 

if  new_lower  >  data^lower 
new^lower  =  data_lower  ; 
end  ; 

if  new_upper  <  data_upper 
new_upper  =  data_upper  ; 
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end  ; 

%  New  limits  must  maintain  a  buffer  distance  from  the  best 
y,  parameter  value  (  Note  -  the  best  value  can  thus  push 
%  the  limits  outward  !  ) 


7o 

% 


if  new_lower  >  (  best_param  -  squeeze„buff er  ) 
new_lower  -  best_param  -  squeeze.buff er  ; 
end  ; 

if  new.upper  <  (  best.param  +  squeeze^buff er  ) 
new^upper  =  best^param  +  squeeze_buffer  ; 
end  ; 

y«  Round  to  integer  value  (  toward  +  /  -  infinity  ) 

new^upper  =  ceil(  new_upper  )  ; 

new^lower  =  floor (  new^lower  )  ; 

7o  Do  not  let  limits  exceed  initial  values 

if  new^lower  <  init_scale(  np,  1  ) 
new^lower  =  init„scale(  np,  1  ); 
end  ; 

if  new.upper  >  init_scale(  np,  2  ) 
new^upper  =  init_scale(  np,  2  ); 
end  ; 

out„put  =  [  igen,  np,  best_param,  new.lower,  new_upper  ] 
pause 

kg_scale(  np,  1  )  *  new^lower  ; 
kg.scaleC  np,  2  )  =  new.upper  ; 

%  Recompute  chromosome  strings  based  on  new  scale  values 


for  ipo  =  1  :  popsize 


new_string  = 

if  np  ==  1 
old_chrom( 
elseif  np  == 
old_chrom( 
elseif  np  == 
old_chrom( 
end  ; 


conv2str(  old_param(  ipo,  np  ) ,  .. 

new^lower,  new^upper 

ipo,  1:16  )  =  new.string  ; 

2 

ipo,  17:32  )  =  new^string  ; 

3 

ipo,  33:48  )  =  new_string  ; 


end  ;  %  for  ipo  loop 


end  ;  %  for  np  loop 
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end  ;  %  End  of  Squeeze  limits 

%  Create  New  Generation 

%  Copy  top  ranked  individuals  exactly  into  next  generation 
if  topsize  >  0 

[  rank_fitnes8,  rank_index  ]  =  sort(  old_fitness  )  ; 

top_chrom  =  old_chrom(  rank_index(  popsize  -  topsize  +  1  . . 

:  popsize  ) ,  :  )  ; 

new_chrom  =  [  new_chrom  ;  top_chrom  ]  ; 
end  :  %  (  if  topsize  ) 

%  Generate  number  of  new  random  individuals  in  next  generation 
if  newsize  >  0 

random_chrom  =  fix(  1.9999999  *  rand(  newsize,  stringlength  )  ); 
new_chrom  =  [  new_clirom  ;  random_chrom  3  ; 
end  ;  %  (  if  newsize  ) 

%  Generate  remaining  number  of  individuals  by  mating  parents 
%  of  previous  generation 

nloops  =  fix(  (  popsize  -  topsize  -  newsize  )  /  2  )  ; 
for  i  =  1  :  nloops 
y,  -  Select  Parents 

[  p_l,  p_2  ]  =  select(  popsize,  sum_fitness,  old.fitness  )  ; 

parent_l  =  old_chrom(  p_l ,  :  )  ; 

parent_2  =  old_chrom(  p_2,  :  )  ; 

y,  -  Create  2  Children,  and  Perform  Crossover 

chlld_l  =  parent_l  ; 

child_2  =  parent_2  ; 

y,  -  Perform  Crossover  separately  on  each  parameter 
for  ip  =  1  :  numparams 
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if  rand  <=  pcross 
if  xover_flag  ==  1 

xsite  =  selxsite(  paramlength,  gen  )  ; 
else 

xsite  =  fix(  rand  *  (  paramlength  -  1  ))  +  1  ; 
end  ; 

xl  =  (  ip  -  1  )  *  paramlength  +  1  ; 
x2  =  xl  +  xsite  -  1  ; 

child.K  xl  :  x2  )  =  parent _2(  xl  :  x2  )  ; 

child_2(  xl  :  x2  )  ==  parent_l(  xl  :  x2  )  ; 

end  ;  %  (  if  rand  ) 

end  ;  %  (  for  ip  ) 

%  Perform  bit  by  bit  mutation  on  each  child 
mu^string  =  rand(  1,  stringlength  )  ; 
index  *  f ind(  mu_string  <«  pmutation  )  ; 
child_l(  index  )  =  *  child.K  index  )  ; 
mu_string  =  rand(  1,  stringlength  )  ; 
index  =  find(  mu.string  <=  pmutation  )  ; 
child_2(  index  )  =  "'  child_2(  index  )  ; 
new^chrom  =  [  new_chrom  ;  child_l  ;  child_2  ]  ; 
end  ;  %  (  For  i  ) 


%  Evaluate  Fitness  of  New  Population 

[new_f itness,  new^param]  =  evalpop7(  A,  B,  C,  D,  ts,  dts,  tl,  dtl, 
hoe.thresh,  phi_rise,  act^thresh,  popsize,  kg_scale,  new^chrom  ) 

%  Option  to  combine  old  population  with  new  population  and  only 
7,  save  best  individuals  from  both 

if  combine .flag  ==  1 

combine.f itness  =  [  old.fitness  ;  new.fitness  ]  ; 
combine.chrom  =  [  old.chrom  ;  new.chrom  ]  ; 
combine.param  =  [  old.param  ;  new.param  ]  ; 
combine.size  =  2  *  popsize  ; 

[  rank.f itness,  rank.index  ]  =  sort(  combine.! itness  )  ; 

new.index  =  rank_index(  combine.size  -  popsize  +  1  :  , . 

combine. size  )  ; 
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new^chrora  =  combiiie_chrom(  [  new_index  ]  ,  :  )  ; 

new^param  =  combine_param(  [  new_index  ] ,  :  )  ; 

new_fitness  =  combine_f itness(  [  new^index  ]  )  ; 

7,  Debug  output  : 

%  rank_f itness,  pause 
%  rank^index,  pause 
7*  combine.fitness,  pause 
7o  new_fitness,  pause 
%  combine ^pararn,  pause 
7o  new_param,  pause 

end  ; 

%  Evaluate  Population  Statistics 
avg^fitness  =  sum(  new.iitness  )  /  popsize  ; 

[  max^fitness,  iraax  ]  "  max(  new_fitness  )  ; 
max^param  =  new_param(  imax,  :  )  ; 

max_chrom  =  new_chrom(  imax,  :  )  ; 

[  min_fitness,  imin  ]  =  min(  new_fitness  )  ; 
min^param  =  new_param(  imin,  :  )  ; 

min^chrom  =  new_chrom(  imin,  :  )  ; 

7o  Save  generational  data  : 

gen^avg.f itness  =  [  gen^avg.f itness  ;  avg^fitness  ]  ; 
gen.max^f itness  =  [  gen_max_f itness  ;  max_fitness  ]  ; 
gen.max^param  =  [  gen_max_param  ;  max^param  ]  ; 

gen_max_chrom  =  [  gen.max.chrom  ;  max_chrom  ]  ; 

gen^min.f itness  =  [  gen.min^f itness  ;  min.fitness  ]  ; 
gen_min.param  =  [  gen.min.param  ;  min.param  ]  ; 

gen.min.chrom  =  [  gen.min.chrom  ;  min.chrom  ]  ; 

gen.kg.lower  =  [  gen.kg.lower  ;  kg_scale(  :,!).*]  ; 
gen.kg.upper  =  [  gen.kg.upper  ;  kg.scale(  :,  2  ).’  ]  ; 
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%  Display  or  plot  values  to  screen 
if  cont_plot  ==  1 
clg  ; 

%  Plot  of  fitness  vs  individual  for  this  generation  : 

subplot (221)  ; 

ind.num  =  1  :  popsize  ; 

axis  ;  %  reset  scaling  to  automatic 

plot(  ind^num,  new_fitness  ), 

title ([’Generation  Number  *,  int2str(gen)] 
xlabel ( ’ Individual ’ ) , . . 
ylabel( ’Fitness’ ) , . . 
grid 


%  Plot  parameter  kl  vs  individual  for  this  generation  ; 

subplot (223)  ; 

ind_num  =  1  :  popsize  ; 

axis([  0,  popsize,  kg_scale(l ,1) ,  kg_scale(l ,2)  ]  )  ; 
plot (  ind^num ,  new_param( : , 1)  ) , 

title ( [’Generation  Number  int2str(gen)]  ),., 
xlabel ( *  Individual ’ ) , . . 
ylabel ( ’ Parameter :  kl ’ ) , . . 
grid 

%  Plot  parameter  k2  vs  individual  for  this  generation  : 

subplot (222)  : 

ind.num  =  1  :  popsize  ; 

axis([  0,  popsize,  kg_scale(2,l) ,  kg_scale(2,2)  ]  )  ; 
plot(  ind_num,  new_param( : ,2)  ), 

t it le( [’Generation  Number  ’,  int2str(gen)]  ),.. 
xlabel ( ’ Individual ’ ) , . . 
ylabel ( ’Parameter :  k2’),.. 
grid 

%  Plot  parameter  k3  vs  individual  for  this  generation  : 

subplot (224)  ; 

ind^num  -  1  :  popsize  ; 

axis([  0,  popsize,  kg_scale(3,l) ,  kg_scale(3,2)  3  )  ; 
plot (  ind_num ,  new.par am ( : , 3)  ) , 

title ( [’Generation  Number  ’ ,  int2str(gen)]  ) , . , 
xlabel ( *  Individual ’ ) , . . 
ylabel ( ’ Parameter :  k3  * ) , . . 
grid 


else 

generation  =  [  gen,  avg^fitness,  max^fitness,  min_fitness  ] 
end  ; 


end  ;  %  (for  igen  loop  ) 
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if  auto_plot  ==  1,  print,  shg,  end  ;  %  output  to  printer 

%  char  =  input (  ’  To  continue  with  more  generations  enter  :  y  ’s’ 

) 

%  if  strcmpC  char,  ’y’  )  ==  1 

%  shg  : 

%  else 

%  break  ; 

7,  end  : 

if  gen  >=  maxgen,  break,  end;  7,  break  operator  control  loop 
end  ;  7*  (  while  -  operator  control  loop  ) 


subplot  7.  return  to  default,  i.e.  full  screen  plots 


%  Plot  of  Maximum  Fitness  Function  : 
gen_num  =  1  :  gen  ; 

axis  ;  %  reset  scaling  to  automatic 

plot(  gen_num,  gen_max_f itness  ), 

title (’Maximum  Fitness  vs.  Generation’),.. 
xlabeK ’Generation’)  . .  . 
ylabelC ’M2iximum  Fitness’ )  , . . 
grid 


%  Always  output  Maximum  Fitness  Function  : 
print  : 


%  Plot  of  best  parameter  vs  generation  number  : 

gen.num  =  1  :  gen  ; 

genpl  =  gen  +  1  ; 

axis([  0.  genpl,  -100,  20  ]  )  ; 
plot(  gen_num,  gen_max_param ( 

gen_num ,  gen_max_param (  :, 2), 
gen_num ,  gen_max_param(  3) . 

title (’Best  Parameters  vs.  Generation’),.. 
xlabeK  ’Generation’) , . . 
ylabel(’Best  Parameters;  kl ,  k2,  k3’),.. 
text(  gen  +  .25  ,  gen_max_param(gen,l) ,  ’kl’), 

text(  gen  +  .25  ,  gen_max_param(gen,2) ,  ’k2’), 

text(  gen  +  .25  ,  gen_max_param(gen,3) ,  ’kS’), 

grid 


7.  Always  output  best  parameters  function  : 


s<»s»«sas« 
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print  : 

%  Output  linal  best  parameters  and  fitness  data  to  screen  : 
gen 

gen_max_f itne8s(  gen  ) 
gen_max_param(  gen,  :  ) 
pause 


%  Plot  of  Average  Fitness  Fiinction  : 

gen_num  =  1  :  gen  ; 

axis  ;  %  reset  scaling  to  automatic 

plot(  gen_num,  gen_avg_f itness  ), 

title (’Genetic  Algorithm  -  Average  Fitness'),.. 
xlabeK  ’Generation’)  , . . 
ylabel (’ Average  Fitness’ ) , . . 
grid, pause 

y.  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’.’s’)  ; 

if  strcmpC  char,  ’y’  )  ==  1 
print 
end  ; 


%  Plot  of  kg  scale  limits  vs  generation  : 
gen_num  =  1  :  gen  ; 

plot (  gen.num ,  gen_kg_lower( : , 1) ,  ’ - ’ ,  gen_num ,  gen_kg_upper ( : , 1)  ,  ’ -  ’ , 
gen_num ,  gen_kg_lower ( : , 2) , ’ - - ’ ,  gen.num ,  gen_kg_upper ( : , 2) . ’ - - ’ , 
gen.num ,  gen.kg_lower( : , 3) , ’ - . ’ ,  gen.num ,  gen.kg.upper ( : , 3)  ,  ’ - .  ’ ) 
title(’-  Gain  Scale  Limits  -’),.. 
xlabeK  ’Generation’ )  , .  . 
ylabel ( ’ Limits ’ ) , . . 
grid,  pause 

y.  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,’s’)  ; 

if  StrcmpC  char,  ’y’  )  ==  1 
print 
end  : 

%  Set  the  gain  to  the  maximum  value  of  the  last  generation  : 

kgain  =  gen.max_param(  gen,  :  ) 

%  Use  program  AP6  to  evaluate  these  gains  ! 
pause 

clc 

%  APGEN3  -  program  completed  ! 
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A.2  AP5 


This  MATLAB  program  builds  the  airframe  and  linear  feedback  controller  state 
space  model  and  provides  several  functions  for  off-line  analysis  and  plots.  It  uses 
the  Control  System  Toolbox  LQR2  fuction  to  determine  the  optimum  steady  state 
Linear  Quadratic  Regulator  design  by  solving  the  associated  continuous-time  Ric- 
cati  equation.  It  is  also  used  to  analyze  other  designs  by  manually  overriding  the 
LQR  gains. 

This  program  then  generates  a  step  response  of  the  resulting  closed  loop 
control  system  and  provides  analysis  of  various  cost  and  fitness  measures.  Next,  it 
generates  step  response  plots  of  the  state  and  control  variables.  Finally,  it  provides 
frequency  domain  analysis  and  plots  of  the  open  loop  control  system. 

A  listing  of  this  program  follows. 


echo  on 
clc 

y,  AP6  -  Autopilot  Design  and  Analysis  Program  #  6 
%  -  airframe  model  only 

%  -  with  linear  feedback  controller 

%  -  for  LQR  design 

y,  by  R.  Hull.  10/6/92 

y. 

%  -  corrected  for  N_delta  sign  error  in  B  matrix 

y.  1/28/93 

%  -  modified  to  use  the  step2  step  response  function 

y,  2/1/93 

%  -  modified  to  use  the  step3  response  function  2/4/93 

y,  -  back  to  step2  function  2/8/93 

%  -  use  step2s  to  compute  dual  dt  step  response  3/1/93 

y, 

y.  This  program  builds  the  linearized  missile  airframe  model 
%  with  unity  accelerometer  and  gyro  feedback  models. 

y, 

y. 

%  Special  Constants  : 

HztR  =  2*pi;  %  Converts  Hz  to  Radians  per  Second 


%  Flight  Condition  Constants  : 
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M. alpha  =  64.11; 

%  M.delta  =  -20.00 

M_delta  =  -62.34; 
M_thetadot  =  0; 

N. alpha  =  .1803; 
N.delta  =  .0738; 
Tau_act  =  .02  ; 
Velocity  =  892.0; 
Flight_Time  =  30.0; 

%  Units  Conversions 


%  Reduced  fin  effectiveness 


%  First  order  actuator  time  constant 
%  meters/sec 
%  seconds 


format  compact  ; 


%  First  -  create  plant  model  without  computation  delay  and  generate 
%  LQR  controller  gains 

%  Define  Airframe  *  First  Order  Actuator  State  Space  Realization  : 

A  =  [  -  N.alpha  1.0  -  N.delta 

M_alpha  0  M_delta 

0  0  -  (  1  /  Tau_act  )  ]  ; 

B  =  [  0 

0 

(  1  /  Tau_act  )  ] 

C  =  [  (  N_alpha  *  Velocity  )  0  (  N_delta  *  Velocity  )  ]  ; 


D  =  C  0  3  ; 


%  Determine  LQR  design  for  basic  plant  with  first  order  actuator  : 

qmatrix  =  [  10  0  0 

0.10 
0  0  .01  ]; 

rmatrix  =  .01; 

%  kgain  =  lqr2 (A , B , qmatrix , rmatrix) ; 

%  LQR  design  gains  : 

%  kgain 
%  pause 


%  kgain  =  [  -34.1723  -3.6403  2.3434  ]  ; 

y.  Genetic  Algorithm  Gains: 

%  kgain  =  [  -48.9113  -6.2638  1.6400  ]  ; 

y,  kgain  =  [  -66.0487  -14.9935  9.1844  ]  ; 

y,  kgain  =  [  -24.7425  -7.0481  1.6289  ]  ; 


%  LQR  gains 

y,  J  Controller 
%  J,Pd  Controller 
%  J.Pd.Psettle  Controller 
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kgain  =  [  -22.3987  -4.5464  1.6514  ]  ;  %  J , Pd, Psettle .Prise  Controller 

Now  -  formulate  open  loop  control  system  lor 
linear  state  feedback  control 

and  determine  gain  and  phase  margins  of  open  loop  system 
Aopen  =  A 
Bopen  =  B 
Copen  =  -  kgain 
Dopen  =  [  0  ] 

Now  -  formulate  closed  loop  regulator  with 
linear  state  feedback  control 

to  determine  step  response  of  closed  loop  system 
hoec  =  1,0  7o  hang  off  error  compensation 

Ac  =  A  -  B  *  kgain 
Be  =  hoec  *  B  *  kgain 
Cc  =  C  -  D  *  kgain 
Dc  =  hoec  *  D  *  kgain 
7o  Set  Cc  k  Dc  matrices  to  force  output  equal  to  control 
Cc  =  -  kgain 
Dc  =  hoec  *  kgain 

y«  Define  hard  limits  for  state  responses  : 

%  (  used  in  step3  function  ) 

Slim  =  [  -  Inf  Inf 

-  Inf  Inf 

-  100.0  100.0  ]  ; 

%  Gains  used  : 
kgain 

pause  %  Strike  any  key  for  step  response  . . . 

%  Now  compute  the  step  response  : 
dtl  =  .01  ; 


tl^final  =5.0  ; 
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dts  =  .001  ; 

ts.final  =0.2  ; 

ts  =  0  :  dts  ;  ts^final  ; 

tl  =  ts^final  :  dtl  :  tl.final  ; 

7o  Now  compute  the  step  response  : 

Rc  =  [  1 

0 

0  ]  ;  %  Reference  Input  Commands 

[ys,yl,xs,xl]  =  8tep2s(Ac ,Bc ,Cc ,Dc ,Rc ,dts,ts,dtl,tl) ; 

%  Evaluate  the  cost  function  separately  for  small  and  large  dt 
%  responses  and  sum  to  get  total  : 

r=[100];  %  indicates  desired  step  response 

jl  =  jcost(  xs,  r,  ys,  qmatrix,  rmatrix,  dts  ); 

j2  =  jcost(  xl,  r,  yl,  qmatrix,  rmatrix,  dtl  ); 

j  =  jl  +  j2  ; 

7o  Evaluate  the  rise  time 

lens  =  length(  ts  )  ; 

t_rise  =  risetime (  xs(:,l),  ts,  .8  )  ; 

if  t_rise  >=  ts(  lens  ) 

t„rise  =  risetime(  xl(:,l),  tl,  .8  )  ; 
end  ; 

7o  Determine  infinity  norm  of  actuator  response 
act^norm  =  infnorm(  xs(:,3)  )  ; 

%  Determine  the  Hang  Off  Error  in  state  1  (  final  value  error  )  : 

[  xrows,  xcols  ]  =  size(  xl  )  ; 

hang_off _error  =  xl(  xrows,  1  )  -  r(  1  )  ; 

%  Gains  used  : 
kgain 

%  Cost  function  : 

j 


7o  Rise  Time  : 
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t_rise 

%  Settling  error  : 
hang_oli_error 

%  Actuator  peak  response 
act_norm 


%  Strike  Enter  to  Actuator  Step  Response  . . . 
pause  ; 

clg 

axis( [0,ts_f inal, -10, 10] ) 

plot(ts,xs(: ,3)) ,  title ( ’Actuator  Step  Response’),.. 
xlabeK ’Time  -  sec’),.. 
ylabelC ’Actuator  Angle  -  deg’),., 
grid, pause 

%  Prompt  lor  bard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,’B’) 

il  strcmpC  char,  ’y’  )  ==  1 
print 
end 


%  Steady  State  Step  Angle  of  Attack  Response  : 
clg 

axis( [0,tl_f inal, 0,2] ) 
plot(ts,xs(: ,1) ,tl,xl(: ,1)) , . . 

title (’Angle  of  Attack  Step  Response’),.. 
xlabelC’Time  -  sec’),.. 
ylabeK ’Angle  of  Attack  -  deg’),., 
grid, pause 

y.  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,’s’) 

if  strcmp(  char,  ’y’  )  ==  1 
print 
end 

%  Strike  any  key  -  for  Gain  and  Phase  Margins  . . . 
pause 

echo  off 

w  =  logspace(-l,7,160) ; 

[mag, phase]  =  bode(  Aopen , Bopen , Copen , Dopen , 1 , w) ; 
[LGm,HGm,Pm,WLg,WHg,Wcp]  =  margin2(  mag,  phase,  w  ) 

y.  Note  -  MARGIN2  is  just  a  modified  form  of  MARGIN  to  accept 
%  mag  in  dB. 
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f printf (  . . 

•  Freq  =  %8.2f  rad/sec  Low  Freq  Gain  Margin  =  yo8.2f  dB\n*,.. 

WLg,  LGm  ) 

f printf (  . . 

’  Freq  =  yo8.2f  rad/sec  High  Freq  Gain  Margin  =  %8.2f  dBXn' , . . 
WHg,  HGm  ) 

f printf (  . . 

’  Freq  =  yo8.2f  rad/sec  Phase  Margin  =  y,8.2f  deg\n’ ,  Wcp,  Pm  ) 
echo  on 


%  Strike  any  key  -  for  Bode  Plot  . . . 
pause 

clg 

subplot (211) 

axis ([-1.7. -100, 100]) 

semilogx(w.mag) ,  title(’Open  Loop  CAS  delta  frequency  response’),.. 
xlabeK  ’Radians  /  sec  ’ )  ,  .  . 
ylabel(’Gain  -  dB’)... 
grid, , . 

subplot (212) 
axis([-l,7,-200,200]) 

semilogx(w, phase) ,  title(*Phase  response’) 

xlabeK  ’Radians  /  sec  * )  , .  . 
ylabel( ’Phase  -  Deg’),., 
grid, . . 

pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,*s’) 

if  strcmp(  char,  ’y’  )  ==  1 
print 
end 


subplot  %  return  to  default,  i.e.  full  screen  plots 

%  Strike  any  key  -  for  Nyquist  plot  . . . 
pause 


w  =  logspace(-l,3,80) ; 

[re.im]  =  nyquist(  Aopen, Bopen, Copen, Dopen, l,w) ; 

clg 

axis 

plot (re, im),  title (’Open  Loop  CAS  delta  Nyquist  Plot’),.. 
xlabeK ’real  axis’),.. 
ylabeK  ’  imaginary  axis  ’ )  ,  .  . 


grid. . . 


pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (•  To  make  hardcopy  of  plot  enter:  y'.’s') 

if  strcmpC  char,  *y’  )  ==  1 
print 
end 


%  Strike  any  key  -  for  Zoomed  Nyquist  plot  : 
pause 


clg 

axis ( [- 10 , 0 , -5 , 5] ) 

plot (re, im),  title (’Open  Loop  CAS  delta  Nyquist  Plot’),. 
xlabelCreal  axis’)  ,  .  . 
ylabel ( ’ imaginary  axis ’ ) , . . 
grid. . . 

pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,’s’) 

if  strcmp(  char,  ’y’  )  ==  1 
print 
end 


clc 

%  APS  -  program  completed  ! 
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A.3  APNMOPT 


This  MATLAB  program  determines  the  true  optimum  fitness  function  for  a  given 
controller  specification,  using  the  MATLAB  function  FMINS,  which  implements  a 
Nelder-Meade  numerical  search  algorithm. 

A  listing  of  this  program  follows. 


echo  on,  clc 
'/.  APNMOPT.  M 
% 

7. 

7. 

% 

7o 

i  by  R.  Hull.  4/6/93 

kg  =  [  -48.9113  -6.2638  1.64  ] 

±it_test  =  -  apfitlunC  kg’  ) 
kg_opt  =  fmins(  ’aplitfun’,  kg’  ) 
lltness  =  -  aplitlun(  kg_opt  ) 
pause  %  APNMOPT  -  program  completed  ! 


-  Autopilot  Nelder-Meade  Optimizaton  Program 

-  with  linear  feedback  controller 

-  numerical  optimization  of  autopilot  controller 

fitness  function  using  Nelder-Meade  simplex 
algorithm  in  Matlab  function  FMINS. 


100 


A.4  APOPLOT 

This  MATLAB  program  plots  the  optimum  fitness  functions  parametrically  vs  one 
gain,  while  holding  the  other  two  constant.  It  overlays  the  genetically  optimized 
fitness  function  with  the  true  optimum  fitness  function  for  comparison. 

A  listing  of  this  program  follows. 


echo  on 
clc 

%  APOPLOT. M  -  Program  to  plot  optimum  fitness  function 

% 

%  by  R.  Hull.  4/6/93 

% 

format  compact  ; 

ga.opt  =  [  -22.3987  -4.6464  1.6614  ]  %  Input  GA  Optimum  Gains 

ga_fit  =  1.2367  %  Input  GA  max  fitness 

true.opt  =  [  -22.2218  -4.4465  1.4819  ]  %  Input  True  Optimum  Gains 

true _f it  =  1.3316  %  Input  True  max  fitness 

y.scale  =  1.8  ;  %  Meuc  y  axis  scale  for  plot 

%  Note  -  change  controller  spec  in  APFITFON  to  match  gains  above  ! 
pause 

kl_low  =  -30  ; 
kl.high  =  -10  ; 
kl_delta  =  1  ; 

k2_low  =  -11  ; 
k2_high  =  -1  : 
k2_delta  =  .6  ; 

k3_low  =  0.1  : 
k3_high  =4.1  ; 
k3_delta  =  .2  ; 

y,  Option  to  skip  kl  loop 

char  =  input (■  To  skip  kl  plot  enter  y;  y’.’s’) 

if  strcmp(  char,  'y*  )  ==  0 
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k2_ga  =  ga_opt (  2  ) ; 
k3_ga  =  ga_opt (  3  ) ; 
k2_true  =  true_opt(  2  ); 
k3_true  =  true_opt(  3  ); 


y,  lor  fixed  value  of  k2  gain 
y,  for  fixed  value  of  k3  gain 

y,  for  fixed  value  of  k2  gain 
y,  for  fixed  value  of  k3  gain 


%  Create  kl  history  to  contain  values  at  ga  optimum  and  true  optimum 
kl_templ  =  kl_low  :  kl_delta  :  kl_high  ; 
kl.insert  =  ga_opt(l)  ; 

kl_temp2  =  [  kl.tempK  kl.templ  <  kl.insert)  kl.insert  .. 
kl_templ(  kl.templ  >  kl.insert)  ]  ; 

kl.insert  =  true.opt(l)  ; 

kl.hist  =  [  kl.temp2(  kl.temp2  <  kl.insert)  kl.insert  . . 
kl_temp2(  kl_temp2  >  kl.insert)  ]  ; 

%  Initialize  for  ki  loop 
max.fitness  =0.0  ; 
ga.fit.hist  =  [  ]  ; 
true.f it.hist  =  [  ]  ; 
kl.len  =  length(  kl.hist  )  : 


for  ki  =  1  :  kl.len 
kl  =  kl.hist (  ki  )  ; 
kg_ga  =  [  kl  k2.ga  k3.ga  ]  ; 
kg.true  =  [  kl  k2.true  k3_true  ]  ; 
ga.fitness  =  -  apfitfun(  kg.ga’  )  ; 
true.f itness  =  -  apfitfun(  kg.true’  )  ; 

%  Save  data  for  plots 

ga.fit.hist  =  [  ga.fit.hist  ;  ga.fitness  ]  ; 

true.fit.hist  =  [  true.f it.hist  ;  true.fitness  ]  ; 

data.out  =  [  kl  ga.fitness  true.fitness  ]  %  output  to  screen 
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end  y,  for  ki  loop 
kl 

y.  Strike  any  key  for  plot  of  Fitness  Function  : 
pause 

clg 

axis([  kl_low,  kl_h.igh,  0,  y_scale  ]  )  ; 

plot(  kl.hist,  ga_fit_hist,  ga_opt(l),  ga_fit.  ’x* .  .. 

kl_hist,  true_lit_hist,  true_opt(l),  true_fit,  '** 
title( ’Fitness  Function  vs  kl'),.. 
xlabeK'kl  gain’),., 
ylabel ( ’ Fitness ’ ) , . . 

text(  ga_opt(l),  ga_fit  +  .1,  ’GA  Optimum’  ),.. 
text(  true.optd),  true _f it  +  .16,  ’True  Optimum’  ),.. 
grid, pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’.’s’) 

if  Btrcmp(  char,  ’y’  )  ==  1 
print 
end 


end  %  end  if  -  option  to  skip  kl  loop 


%  Option  to  skip  k2  loop 
char  =  input ('  To  skip  k2  plot  enter  y:  y’,’s’) 
if  strcmpC  char,  ’y’  )  ==  0 

kl_ga  =  ga_opt(  1  );  %  for  fixed  value  of  kl  gain 

h3_ga  =  ga_opt(  3  );  %  for  fixed  value  of  k3  gain 

kl_true  =  true_opt(  1  );  %  for  fixed  value  of  kl  gain 

k3_true  =  true_opt(  3  );  %  for  fixed  value  of  k3  gain 

y.  Create  k2  history  to  contain  values  at  ga  optimum  and  true  optimum 
k2_templ  =  k2_low  :  k2_delta  :  k2_high  ; 
k2_insert  =  ga_opt(2)  ; 

k2_temp2  =  [  k2_templ(  k2_templ  <  k2_in8ert)  k2_insert  .. 
k2_templ(  k2_templ  >  k2_insert)  ]  ; 


k2_insert  =  true_opt(2)  ; 
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k2_hist  =  [  k2_temp2(  k2_temp2  <  k2_iiisert)  k2_insert  . . 
k2_temp2(  k2_temp2  >  k2_insert)  ]  ; 

%  Initialize  for  ki  loop 
max_fitness  =  0.0  ; 
ga_lit_hLi8t  =  [  ]  : 
true_iit_hist  =  [  ]  ; 
k2_len  =  length (  k2_hist  )  ; 


for  ki  =  1  :  k2_len 
k2  =  k2_hist(  ki  )  ; 
kg-ga  =  C  kl_ga  k2  k3_ga  ]  ; 
kg_true  =  [  kl.true  k2  k3_tnie  ]  ; 
ga_fitness  =  -  apfitfun(  kg_ga’  )  ; 
true .fitness  =  -  apfitfunC  kg.true’  )  ; 

%  Save  data  for  plots 

ga_fit_hist  =  [  ga_fit_hist  ;  ga.fitness  ]  ; 
tme_fit_hist  =  [  tme_fit_hlst  ;  true.fitness  ]  ; 
data.out  =  [  k2  ga.fitness  true.fitness  ]  %  output  to  screen 

end  %  for  ki  loop 
k2 

y,  strike  any  key  for  plot  of  Fitness  Function  : 
pause 

clg 

axis([  k2_low,  k2_high,  0,  y.scale  ]  )  ; 
plot(  k2_hist,  ga_fit_hist,  ga_opt(2),  ga.fit,  ’x’ ,  .. 
k2_hist.  true_fit_hist,  true_opt(2).  true.fit, 
title( 'Fitness  Function  vs  k2’),.. 
xlabel ( ' k2  gain ’ ) , . . 
ylabel ( ’ Fitness ’ ) , . . 

text(  ga_opt(2),  ga.fit  +  .1,  'GA  Optimum'  )... 
text(  true_opt(2),  true.fit  +  .15.  'True  Optimum’  ),.. 
grid, pause 
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%  Prompt  for  hard  copy  plot; 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y'.’s’) 

if  strcmp(  char,  'y'  )  ==  1 
print 
end 


end  %  end  if  -  option  to  skip  k2  loop 


%  Option  to  skip  k3  loop 
char  =  input ('  To  skip  k3  plot  enter  y:  y’.'s’) 
if  strcmp(  char,  ’y’  )  =*  0 


kl_ga  =  ga_opt(  1  ); 
k2_ga  =  ga_opt(  2  ); 
kl_true  =  true_opt(  1  ); 
k2_true  =  true_opt (  2  ) : 


%  for  fixed  value  of  kl  gain 
%  for  fixed  value  of  k2  gain 

*/,  for  fixed  value  of  kl  gain 
%  for  fixed  value  of  k2  gain 


%  Create  k3  history  to  contain  values  at  ga  optimum  and  true  optimum 
k3_templ  =  k3_low  :  k3_delta  :  k3_high  ; 
k3_insert  =  ga_opt(3)  ; 

k3_temp2  =  [  k3_templ(  k3_templ  <  k3_insert)  k3_insert  .. 
k3_templ(  k3_templ  >  k3_insert)  ]  ; 

k3_insert  =  true_opt(3)  ; 

k3_hist  =  [  k3_temp2(  k3_temp2  <  k3_insert)  k3_insert  . . 
k3_temp2(  k3_temp2  >  k3_insert)  ]  ; 


%  Initialize  lor  ki  loop 
max_fitness  =0.0  ; 
ga_lit_hi8t  =  [  ]  : 
true_lit_hist  =  [  ]  : 
k3_len  =  length(  k3_hist  )  ; 


lor  ki  =  1  :  k3_len 
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k3  =  k3_hist(  ki  )  ; 
kg-6^  “  [  kl_ga  k2_ga  k3  ]  ; 
kg_true  =  [  kl_true  k2_true  k3  ]  ; 
ga_litnes8  =  -  apfitlunC  kg_ga’  )  ; 
true  .fitness  =  -  apfitfimC  kg.tnie’  )  ; 

%  Save  data  for  plots 

ga_fit_hist  =  [  ga_fit_hist  ;  ga.fitness  ]  ; 
true_fit_hist  =  [  tnie_lit_hist  ;  true.litness  ]  ; 
data.out  *  [  k3  ga.litness  true.fitness  ]  %  output  to  screen 

end  %  for  ki  loop 
k3 

%  Strike  any  key  for  plot  of  Fitness  Function  : 
pause 

clg 

axis([  k3_low,  kS.high,  0,  y.scale  ]  )  ; 
plot(  k3_hist.  ga_fit_hist,  ga_opt(3).  ga.fit,  ’x’.  .. 
k3_hi8t,  true_fit_hist,  true_opt(3),  true.fit, 
title ( ’Fitness  Function  vs  k3’),.. 
xlabel ( ’ k3  gain ’ ) , . . 
ylabel(’ Fitness’) , . . 

text(  ga_opt(3),  ga.fit  +  .1,  ’GA  Optimum’ 
text(  true_opt(3).  true.fit  +  .15,  ’True  Optimum’ 
grid, pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’.’s’) 

if  strcmpC  char,  ’y’  )  ”  1 
print 
end 


end  %  end  if  -  option  to  skip  k3  loop 
%  APOPLOT  -  program  completed  ! 
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A.5  GATEST 

This  MATLAB  program  implements  a  simple  genetic  algorithm  to  demonstrate  its 
effectiveness  on  a  single  parameter  example  with  a  non-trivial  solution. 

A  listing  of  this  program  follows. 


echo  on 
clc 

%  GATEST  -  Genetic  Algorithm  Test  Program 

%  -  optimize  special  liuiction  of  one  parameter 

% 

%  by  R.  Hull.  2/19/93 

% 

% 

%  Special  Constants  : 

HztR  =  2*pi:  %  Converts  Hz  to  Radians  per  Second 


%  Global  Variables  : 
global  cont_plot  ; 

cont_plot  =  1  :  %  Continuous  plot  flag  (  on  monitor  ) 

auto_plot  =  1  :  %  Automatic  plot  flag  (  to  printer  ) 

combine_flag  =  0  ;  %  Option  to  combine  old  and  new  populations 
%  and  save  the  best  individuals  of  both. 

%  Genetic  Algorithm  Structures 

maxgen  =  4  ;  %  number  of  generations  per  loop 

popsize  =  20  :  %  number  of  individules  in  population 

%  must  be  even  ! 

topsize  =  2  ;  %  number  (  even  )  of  top  ranked  individuals  to 

%  propagate  exactly  in  the  next  generation 

if  combine_flag  ==  1  %  Make  sure  topsize  is  0  if  combine  option 
topsize  =  0  :  %  is  used, 

end  ; 

newsize  =  0  ;  %  number  (  even  )  of  new  random  individuals  to 

%  include  in  each  new  generation 


nuraparams  =  1  ;  %  number  of  parameters  to  optimize 


paramlength  =  16  ;  %  number  of  bits  per  parameter 

stringlength  =  numparams  *  paramlength  ;  %  string  length 

%  (  bits  ) 

t_final  =  5.0  ;  %  fitness  function  final  time  -  sec. 

pcross  =  .8  ;  %  probability  of  crossover 

pmutation  =  .026  %  probability  of  mutation 

kg_scale  =[05  ]  '/,  gain  parameter  scales 

%  Units  Conversions 

%  Initialize  time  values 
dt  =  .01  ; 
t  =  0:dt:t_final  : 

%  Set  screen  display  format  to  suppress  excessive  line  feeds 
format  compact  ; 

%  Initialize  New  Population  to  Random  Strings 

new.chrom  =  fix(  1.9999999  *  rand(  popsize,  stringlength  )  ) 

%  Clear  generational  statistics 
gen_avg_f itness  =  [  ]  ; 

gen_max_f itness  =  [  ]  : 
gen_max_param  =  [  ]  : 
gen_max_chrom  =  [  ]  ; 
gen_min_f itness  =  [  ]  : 
gen_min_param  =  [  ]  : 
gen_min_chrom  =  [  ]  ; 
gen_kg_lower  =  [  ]  : 
gen_kg_upper  =  [  ]  ; 
gen  =  1  ;  %  generation  number 
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y.  Evaluate  Fitness  of  Initial  Population 

[new_fitness,  new.param]  =  evalpoptC  popsize,  kg_scale.  new.chrom  ); 
y.  Evaluate  Initial  Population  Statistics 
avg_litness  =  8um(  new_litness  )  /  popsize  ; 

[  max_fitness,  imax  ]  =  maxC  new_fitness  )  ; 
max.param  =  new_parain(  imax  )  ; 
max_chrom  =  new_chrom(  imax  )  ; 

[  min.fitness,  imin  ]  =  min(  new.fitness  )  ; 
min.param  =  new_param(  imin  )  ; 
min.chrom  =  new_chrom(  imin  )  ; 

%  Save  initial  generation  data  : 

gen_avg_fitness  =  [  gen_avg_f itness  ;  avg.fitness  ]  ; 
gen_max_f itness  =  [  gen_max_f itness  ;  max.fitness  ]  ; 
gen_max_param  =  [  gen_max_param  ;  max.param  ]  ; 

gen.max.chrom  =  [  gen_max_chrom  ;  max.chrom  ]  ; 

gen_min_f itness  =  [  gen_min_f itness  ;  min.fitness  ]  ; 
gen.min.param  =  [  gen.min.param  ;  min.param  ]  ; 

gen.min.chrom  =  [  gen.min.chrom  ;  min.ch.rom  ]  ; 

%  Plot  of  Fitness  Function  vs  x  : 
fit.t  =  [  ]  : 
t  =  0: .01:5.0  ; 
for  i  =  1  :  length (t) 

fit  =  -  (  t(i)  -  2.5  )  “  2  +  10.0  +  5  *  sin(  4  *  pi  *  t(i)  ); 

fit.t  =  [  fit.t  :  fit  ]  : 
end; 

plot(  t.  fit.t  ). 

title ( ’Fitness  Function’),.. 
xlabeK  ’x’ )  , . . 
ylabeK  ’Fitness  ’ )  , .  . 
grid, pause 


%  Prompt  for  hard  copy  plot; 


char  =  input (’  To  make  hardcopy  oT  plot  enter;  y'.’s’) 


if  strcmp(  char,  ’y’  )  ==  1 
print 
end  ; 


%  Display  or  plot  values  to  screen 
pause 

if  cont_plot  ==  1 
clg  : 

%  Plot  of  fitness  vs  individual  for  this  generation  : 

subplot (221)  ; 

ind_num  =  1  :  popsize  ; 

plot(  ind_num,  new_fitness  ), 

title([ ’Generation  Number  int2str(gen)]  ),. 
xlabeK  ’  Individual  ’ )  ,  .  . 
ylabelC ’Fitness’) . 
grid 


%  Plot  of  Maximum  fitness  vs  generation  ; 

subplot (222)  ; 

gen_num  =  1  :  gen  ; 

plot(  gen_num,  gen_max_f itness  ). 

title(’-  Maximum  Fitness  -’),.. 
xlabeK ’Generation’)  , . . 
ylabel( ’Maximum  Fitness’),., 
grid 


%  Plot  parameter  vs  individual  for  this  generation  : 

subplot (223)  : 

ind_num  =  1  :  popsize  ; 

plot (  ind_num ,  new_param  ) . 

title(C’ Generation  Number  int2str(gen)]  ), 
xlabel( ’ Individual ’ ) . . . 
ylabel ( ’ Parameter ;  x ’ ) , . . 
grid 


%  Plot  of  Average  Fitness  vs  generation  : 

subplot (224)  ; 

gen_num  =  1  :  gen  ; 

plot(  gen_num,  gen_avg_f itness  ), 

title(’-  Average  Fitness  -’),.. 
xlabeK  ’Generation’ )  , .  . 
ylabel ( ’ Average  Fitness ’ ) . . . 
grid, pause 
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%  if  auto_plot  ==  1 ,  print ,  shg ,  end  ; 
y,  Prompt  for  hard  copy  plot: 

char  =  input (•  To  make  hardcopy  of  plot  enter:  y’.’s’)  ; 

if  strcmp(  char,  'y’  )  ==  1  .  print,  end  ; 

else 

generation  =  [  gen.  avg_fitness.  max_fitness.  min_fitness  ] 
end  ; 

while  1  %  operator  control  loop 

%  Genetic  Algorithm  Optimization  Loop 
for  igen  =  1  :  maxgen 
gen  =  gen  +  1  : 
y.  Save  Old  Generation 
old_fitness  =  new_fitness  ; 
old_chrom  =  new_chrom  ; 
old_param  =  new_param  ; 
new_chrom  =  [  ]  ; 

8um_fitness  =  sum(  old_fitness  )  ; 
y.  Create  New  Generation 

%  Copy  top  ranked  individuals  exactly  into  next  generation 
if  topsize  >  0 

[  rank_f itness.  rank_index  ]  =  sort(  old_fitness  )  ; 

top_chrom  =  old_chrom(  rank_ index (  popsize  -  topsize  +  1  .  . 

:  popsize  ) .  :  )  ; 

new_chrom  =  [  new_chrom  ;  top_chrom  ]  ; 
end  ;  %  (  if  topsize  ) 

%  Generate  number  of  new  random  individuals  in  next  generation 
if  newsize  >  0 

random_chrom  =  fix(  1 .9999999*rand (newsize ,  stringlength) ) ; 
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new_chrom  =  [  new_chroin  ;  random_chrom  ]  ; 
end  ;  %  (  if  newsize  ) 

%  Generate  remaining  number  of  individuals  by  mating  parents 
%  of  previous  generation 

nloops  =  fix(  (  popsize  -  topsize  -  newsize  )  /  2  )  ; 
for  i  =  1  ;  nloops 
%  -  Select  Parents 

C  P_l.  p_2  ]  =  select (  popsize.  sum.fitness,  old.fitness  )  ; 
parent_l  =  old_chrom(  p_l,  :  )  ; 
parent_2  =  old_chrom(  p_2,  :  )  ; 

%  “  Create  2  Children,  and  Perform  Crossover 

child_l  =  parent_l  ; 
child_2  =  parent_2  ; 

%  -  Perform  Crossover  separately  on  each  parameter 

for  ip  =  1  :  numparams 
if  rand  <=  pcross 

xsite  =  fix(  rand  *  (  paramlength  -  1  ))  +1  ; 
xl  =  (  ip  -  1  )  *  paramlength  +  1  ; 
x2  =  xl  +  xsite  -  1  ; 

child_l(  xl  :  x2  )  =  parent_2(  xl  :  x2  )  ; 

child_2(  xl  ;  x2  )  =  parent_l(  xl  :  x2  )  ; 

end  ;  %  (  if  rand  ) 

end  :  %  (  for  ip  ) 

y,  Perform  bit  by  bit  mutation  on  each  child 
mu_Btring  =  reoidC  1,  stringlength  )  ; 
index  =  find(  mu_string  <=  pmutation  )  ; 
child_l(  index  )  =  "  child_l(  index  )  ; 
mu_string  =  rand(  1,  stringlength  )  ; 
index  =  find(  mu_string  <=  pmutation  )  ; 
child_2(  index  )  =  "  child_2(  index  )  ; 
new_chrom  =  [  new.chrom  ;  child_l  ;  child_2  ]  ; 


end  ;  %  (  For  i  ) 
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%  Evaluate  Fitness  of  New  Population 

[  new_fitness,  new_pararo  ]  *  evalpopt(  popsize,  kg_scale,  .. 

new_chrom  )  ; 

%  Option  to  combine  old  population  with  new  population  and  only 
%  save  best  individuals  from  both 

if  combine_flag  ==  1 

combine .fitness  =  [  old.fitness  ;  new.fitness  ]  ; 
combine.chrom  =  [  old.chrom  ;  new.chrom  ]  ; 
combine _param  =  [  old.param  ;  new.param  ]  ; 
combine.size  =  2  *  popsize  ; 

[  rank.f itness,  rank.index  ]  =  sort(  combine .fitness  )  ; 

new.index  =  rank_index(  combine.size  -  popsize  +  1  :  . . 

combine.size  )  ; 

new.chrom  =  combine.chrom(  [  new.index  ] .  :  )  ; 
new.param  =  combine_param(  [  new.index  ] ,  :  )  ; 
new.fitness  =  combine.fitnessC  [  new.index  ]  )  : 
end  ; 

%  Evaluate  Population  Statistics 
avg.fitness  =  sum(  new.fitness  )  /  popsize  ; 

[  max.fitness,  imax  ]  =  max(  new.fitness  )  ; 
max.param  =  new.param (  imax  )  ; 
max.chrom  =  new.chrom (  imax  )  ; 

[  min.fitness,  imin  ]  =  min(  new.fitness  )  ; 
min.param  =  new.param!  imin  )  ; 
min.chrom  =  new.chrom!  imin  )  ; 

%  Save  generational  data  : 

gen.avg.fitness  =  [  gen.avg.f itness  ;  avg.fitness  ]  ; 


gen_max_litness 
gen_max_param  = 
gen_max_chrom  = 
gen_inin_f  itness 
gen_min_param  = 
gen_min_chrom  = 
geii_kg_lower  = 
gen_kg_upper  = 


=  [  gen_max_f itness  ;  max_f itness 
[  gen_max_param  ;  max_param  ]  ; 

[  gen_max_chrom  ;  max_chrom  ]  ; 

=  [  gen_min_f itness  ;  min_litness 
[  gen_inin_parajn  ;  min_param  ]  ; 

[  gen_min_ch.rom  ;  min_chrom  ]  ; 

[  gen_kg_lower  ;  kg_scale(  1  )  ]  ; 
[  gen_kg_upper  ;  kg_scale(  2  )  ]  ; 


] 


] 


%  Display  or  plot  values  to  screen 
if  cont_plot  ==  1 
clg  ; 

%  Plot  of  fitness  vs  individual  for  this  generation  : 

subplot (221)  ; 

ind_nuin  =  1  :  popsize  ; 

plot(  ind_num,  new_f itness  ), 

title([’ Generation  Number  int2str(gen)]  ),. 
xlabel ( ’ Individual ' ) . . • 
ylabel( ’Fitness’) , . . 
grid 


y.  Plot  of  Maximum  fitness  vs  generation  ; 

subplot (222)  ; 

gen_num  =  1  :  gen  ; 

plot(  gen_num,  gen_max_f itness  ), 

title ( ’ -  Maximum  Fitness  -’),.. 
xlabeK’ Generation’  )  , . . 
ylabel (’ Maximum  Fitness’) , . . 
grid 


%  Plot  parameter  vs  individual  for  this  generation  ; 

subplot (223)  ; 

ind_num  =  1  :  popsize  ; 

plot (  ind_num ,  new_param  ) , 

title([’ Generation  Number  ’.  int2str(gen)]  ). 
xlabeK  ’  Individual  ’ )  , .  . 
ylabel ( ’ Parameter :  x ’ ) , . . 
grid 


y.  Plot  of  Average  Fitness  vs  generation  : 
subplot (224)  : 
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gen_num  =  1  :  gen  ; 
plot(  gen.num,  gen_avg_f itness  ), 
titleC’-  Average  Fitness 
xlabel( 'Generation’) . . . 
ylabelC 'Average  Fitness’) , . . 
grid 


pause  ; 

y,  if  auto_plot  ==  1,  print,  shg,  end  ; 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’.’s’)  ; 

if  strcmp(  char,  'y'  )  ==  1  .  print,  end  ; 


else 

generation  *  [  gen,  avg_fitness,  max.fitness,  min.fitness  ] 
end  ; 

end  :  y,  (for  igen  loop  ) 

gen 

gen_max_f itness (  gen  ) 
gen_max_param(  gen,  :  ) 

char  =  input (’To  continue  with  more  generations  enter  :  y  ’,  ’s’) 

if  strcmp(  char,  ’y’  )  ==  1 
shg  : 
else 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’.’s’)  ; 

if  strcmp(  char,  ’y’  )  ==  1  .  print,  end  ; 
break  ; 
end  ; 

end  ;  %  (  while  -  operator  control  loop  ) 

%  Output  final  best  parameters  and  fitness  data  to  screen  : 
gen 

gen_max_f itness (  gen  ) 
gen_max_param(  gen,  :  ) 
pause 

subplot  %  return  to  default,  i.e.  full  screen  plots 
%  Plot  of  Average  Fitness  Function  : 


gen.num  =  1  :  gen  ; 

plot(  gen.num,  gen_avg_f itness  ), 


title ( ’Genetic  Algorithm  -  Average  Fitness’), 
xlabeK ’Generation’)  , . . 
ylabel( ’Average  Fitness’) , . . 
grid, pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (’  To  make  hardcopy  of  plot  enter:  y’,'s’) 

if  strcmp(  char,  ’y’  )  ==  1 
print 
end  ; 


%  Plot  of  Maximum  Fitness  Function  : 

gen_num  =  1  :  gen  ; 

plot(  gen_num,  gen_max_f itness  ), 

title (’Genetic  Algorithm  “  Maximum  Fitness’), 
xlabeK  *  Generation  ’ )  , . . 
y label ( ’Maximum  Fitness ’ ) , . . 
grid, pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (•  To  make  hardcopy  of  plot  enter:  y’.’s’) 

if  strcmpC  char,  *y*  )  ==  1 
print 
end  ; 


%  Plot  of  best  parameter  vs  generation  number  : 

gen^num  =  1  :  gen  ; 

plot (  gen^num ,  gen_max_param  ) , 

title (’Genetic  Algorithm  -  Best  Parameter’),. 
xlabeK ’Generation’)  , .  , 
ylabeK’Best  Parameter:  x’)  ,  .  . 
grid, pause 

%  Prompt  for  hard  copy  plot: 

char  =  input (*  To  make  hardcopy  of  plot  enter:  y*,*s*) 

if  strcmp(  char,  ’y*  )  ==  1 
print 
end  ; 


clc 

%  GATEST  -  program  completed  ! 
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APPENDIX  B 


Special  MATLAB  Supporting  Functions 


Several  special  supporting  functions  were  written  in  MATLAB  and  are  required  by 
one  or  more  of  the  MATLAB  programs  described  in  Appendix  A.  This  appendix 
briefly  describes  these  supporting  functions  and  provides  a  listing  of  their  MATLAB 
source  files.  In  addition  these  functions  and  the  main  programs  require  numerous 
standard  MATLAB  functions,  as  well  as  supporting  functions  from  the  MATLAB 
Control  System  Toolbox  [13,10]. 

B.l  APFITFUN 

This  MATLAB  function  evaluates  the  fitness  function  for  a  given  set  of  gain  pa¬ 
rameters  and  design  specification.  It  first  calculates  the  system  step  response,  using 
STEP2S,  and  then  determines  the  linear  quadratic  performance  index  using  JCOST. 
It  then  evaluates  the  cost  function  penalties  imposed  by  peak  actuator  response, 
settling  error,  and  rise  time  constraints.  The  final  cost  function  is  limited  and 
converted  to  a  fitness  function  appropriate  for  the  genetic  algorithm. 

A  program  listing  of  this  function  follows. 
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% 

% 

% 

% 

I 
% 

% 

% 

% 

% 

function  fitness  =  apfitfun(  kg_vector  )  ; 


kg  =  kg_vector’  ; 

%  convert  from  column  vector  required  by 
%  function  fmins 

Flight  Condition 

Constants  : 

M_alpha  =  64.11; 
M_delta  =  -62.34; 

M^thetadot  =  0; 

N_alpha  =  .1803; 
N_delta  =  .0738; 

Tau_act  =  .02  ; 

%  Actuator  first  order  time  constant 

Velocity  =  892.0 

;  y,  meters/sec 

%  Performance  Specifications 

%  ph,i_rise  =  6.0  ;  %  -  removes  rise  time  penalty 

phi_rise  =  .30  ;  %  rise  time  specification 

y.  act_thresh  =  1000.0  ;  %  -  removes  actuator  response  penalty 

act_thresh  =  6.0  ;  %  maximum  actuator  response  specification 

%  penalty  threshold 

7.  hoe_thresh  =  100.0  ;  %  -  removes  settling  error  penalty 

hoe.thresh  =  0.0  ;  %  settling  error  penalty  threshold 

%  specification 


%  Define  time  histories  : 
dtl  =  .01  ; 
tl_final  =6.0  ; 
dts  =  .001  ; 
ts_final  =  0.2  ; 
ts  =  0  :  dts  :  ts_final  ; 
tl  =  ts_final  ;  dtl  ;  tl_linal  ; 


FUNCTION  APFITFUN  -  Autopilot  Fitness  Function 

Note:  this  function  duplicates  much 
of  the  code  in  EVALP0P7 
Computes  the  Fitness  Function  based 
on  three  input  gain  parameters  (  kg  ) 

-  Inverts  final  fitness  value  for  use 

by  MATLAB  fmins  optimization  function. 

by  :  R.  Hull  4/6/93 


%  Define  Airframe  +  First  Order  Actuator  State  Space  Realization 


A  =  [ 


N_alpha  1 . 0 
M_alpha  0 
0  0 


-  N_delta 
M.delta 

-  (  1  /  Tau_act  )  ]  ; 


B  =  [  0 

0 

(  1  /  Tau_act  )  ]  ; 

C  =  [  (  N.alpha  *  Velocity  )  0  (  N.delta  *  Velocity  )  ]  ; 


D  =  [  0  3  : 


%  Determine  LQR  design  for  basic  plant  with  first  order  actuator 

qmatrix  =  [  10  0  0 

0.10 
0  0  .01  ]; 

rmatrix  =  .01; 


%  Formulate  closed  loop  regulator  with 
%  linear  state  feedback  control 

to  determine  step  response  of  closed  loop  system 
Set  Cc  ft  Dc  matrices  to  force  output  equal  to  control  for  step 
response 

hoec  =  1.0  ;  %  hang  off  error  compensation 

Ac  =  A  -  B  *  kg  : 

Be  =  hoec  *  B  *  kg  ; 

Cc  =  -  kg  : 

Dc  =  kg  ; 

%  Determine  the  closed  loop  step  response  : 

Rc  =  [  1 

0 

0  3  : 

%  Small  and  large  dt  step  responses  : 

[ys.yl.xs,xl3  =  step2s(Ac ,Bc .Cc .Dc .Rc .dts .ts.dtl.tl) ; 

%  Evaluate  the  cost  function  separately  for  small  and  large  dt 
%  responses  and  sum  to  get  total  : 
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r  =  [  1  0  0  ]  ;  %  indicates  desired  step  response 

jl  =  jcost(  xs,  r,  ys,  qmatrix,  rmatrix,  dts  ); 

j2  =  jcost(  xl,  r,  yl.  qmatrix,  rmatrix,  dtl  ); 

j  =  jl  +  j2  : 

%  Evaluate  the  rise  time 

lens  =  length (  ts  )  ; 

t_rise  =  risetime(  xs(:,l),  ts,  .8  )  ; 

if  t_rise  >=  ts(  lens  ) 

t_rise  =  risetimeC  xl(:,l),  tl,  .8  )  ; 
end  ; 

%  Determine  infinity  norm  of  actuator  response 
act_norm  =  inf norm (  xs(:.3)  )  ; 

%  Determine  the  Hang  Off  Error  in  state  1  (  settling  error  )  : 

[  xrows,  xcols  ]  =  size(  xl  )  ; 

hang_off_error  =  xl(  xrows,  1  )  -  r(  1  )  ; 

if  isnan(  hang_off _error  ) 
hang_off _error  =  10  ; 
end  ; 

if  abs(  hang_off .error  )  >  10 
hang.off .error  =  10  ; 
end  ; 

y.  Convert  cost  function  to  suitably  limited  fitness  function 
%  including  penality  for  exceeding  rise  time  specification 

if  isnan(  j  ) 
j  =  100  ; 
end  ; 

if  j  >  100 
j  =  100  ; 
end  ; 

if  t.rise  >  phi.rise 

cost.term  =  j  +  10.0  +  (  t.rise  -  phi.rise  )  ; 

else 

cost.term  =  j  : 
end  ; 


%  Penalty  for  Settling  Error 

if  abs(  hang_off_ error  )  >  hoe.thresh 

cost.term  =  cost^term  +  10.0  *  abs(  hang_of f _error  ) 

end  ; 

%  Penalty  for  infinity  norm  of  actuator  response 

if  abs(  act_norm  )  >  act^thresh 
cost_term  =  cost_term  . . 

+  20.0  *  abs(  act_norm  “  act_tliresh  ) 

end  ; 

%  Limit  Cost  Term 

if  cost.term  >  100 
cost^term  =  100  ; 

elseif  cost .term  <  .6 
cost.term  =  .5  ; 

end  ; 

fitnessl  =  exp(  2  /  cost.term  )  -  1  ; 

fitness  =  -  fitnessl  ;  %  Invert  for  use  in  FMINS 


%  End  APFITFUN 
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B.2  CONV2NM2 


This  MATLAB  function  converts  a  16  bit  character  string  to  an  appropriately  scaled 
real  number  ranging  between  the  specified  lower  and  upper  limits. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  C0NV2NM2 

%  This  function  converts  the  input  string  variable  (  actually 
%  an  16  column  matrix  of  zero’s  and  one’s  )  to  a  number  (  num  ) 
%  with  value  between  lower  and  upper. 

%  Note  that  in  some  cases  there  may  not  be  an  exact  zero  value 
%  after  the  conversion. 

function  [num]  =  conv2nm2(  string,  lower,  upper  ) 

d  =  [  32768 
16384 
8192 
4096 
2048 
1024 
512 
256 
128 
64 
32 
16 
8 
4 
2 

1  ]  ; 

num  =  lower  +  (  upper  -  lower  )  *  (  string  *  d  )  /  65535  ; 

%  End  of  function  C0NV2NM2 
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B.3  CONV2STR 


This  MATLAB  function  converts  a  real  number  to  a  16  bit  character  string,  assum¬ 
ing  it  ranges  between  the  specified  lower  and  upper  limits. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  C0NV2STR 

%  This  function  converts  the  input  number  to  a  16  bit  string 
%  (  actually  a  16  column  matrix  of  zero’s  and  one’s  ). 

%  It  is  assummed  the  input  number  lies  between  the  lower  and 
%  upper  limit.  Usually  the  input  number  has  been  previously 
%  generated  from  a  string  and  scaled  by  using  conv2nm2. 

function  [string]  =  conv2str(  num,  lower,  upper  ) 

rem  =  66535  *  (  num  -  lower  )  /  (  upper  -  lower  )  ; 

div  =  32768  ; 

string  =  zeros (  1,  16  )  ; 

for  i  =  1:16 

bit  =  fix(  rem  /  div  )  ; 

string (  i  )  =  bit  ; 

rem  =  rem  -  bit  *  div  ; 

div  =  div  /  2  ; 

end  ; 

%  End  of  function  C0NV2STR 
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B.4  EVALPOP7 

This  MATLAB  function  evaluates  the  genetic  algorithm  fitness  function  for  a  given 
system  configuration  and  design  specification.  It  converts  the  parameter  string 
information  into  appropriately  scaled  gain  values,  and  evaluates  the  fitness  function 
for  every  member  of  the  population.  It  first  calculates  the  system  step  response, 
using  STEP2S,  and  then  determines  the  linear  quadratic  performance  index  using 
JCOST.  It  then  evaluates  the  cost  function  penalties  imposed  by  peak  actuator 
response,  settling  error,  and  rise  time  constraints.  The  final  cost  function  is  limited 
and  converted  to  a  fitness  function  appropriate  for  the  genetic  algorithm. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  EVALP0P7  -  Evaluate  Population  function  M  file 
%  -  assuming  three  parameters  per  string  ;  kl , k2 , k3 

%  -  using  LQR  cost  function 

%  -  includes  penality  for  rise  time  specification 

%  -  use  step2  function 

%  -  include  actuator  model 

%  -  include  hang  off  error  penalty 

%  -  include  actuator  max  response  penalty 

%  -  use  smaller  dt  to  calculate  actuator  response  norm 

i  by  R.  Hull.  3/1/93 

% 

% 

function  [  pop_fitness,  pop_param  ]  =  evalpop7(  A,  B,  C,  D,  ts,  dts,  .. 
tl,  dtl,  hoe_thresh.  phi_rise ,  act_thresh,  popsize ,  kg_scale,  .. 
pop_chrom  ) ; 

pop_param  =  [  ]  ; 

pop_fitness  =  [  ]  ; 

kg  =  [  1.0  1.0  1.0  ]  : 

%  Set  cost  function  weighting  matrices 

qmatrix  =  [  10  0  0 

0  .1  0 
0  0  .01  ]; 


rmatrix  =  . 01 ; 


%  Population  loop 

for  k  =  l:popsize 

kg(l)  =  conv2nm2(  pop.chrom (k ,  1:16),  kg_scale(l , 1) , . . 

kg_scale(l ,2)  )  ; 

kg(2)  =  conv2nin2(  pop_chrom(k,  17:32),  kg_scale(2 ,1)  ,  .  . 

kg_scale(2,2)  )  ; 

kg(3)  =  conv2nm2(  pop_chrom(k,  33:48),  kg_scale(3,l) , . . 

kg_scale(3,2)  )  ; 

pop_param  =  [  pop_parani  ;  kg  ]  ; 

%  Formulate  closed  loop  regulator  with 
linear  state  feedback  control 

to  determine  step  response  of  closed  loop  system 
Set  Cc  k  Dc  matrices  to  force  output  equal  to  control 

hoec  =  1.0  :  %  hang  off  error  compensation 

Ac  =  A  -  B  *  kg  ; 

Be  =  hoec  *  B  *  kg  : 

Cc  =  -  kg  : 

Dc  =  kg  ; 

%  Determine  the  closed  loop  step  response  : 

Rc  =  [  1 

0 

0  ]  : 

%  Small  and  large  dt  step  responses  : 

[ys,yl,xs,xl]  =  step2s(Ac ,Bc,Cc ,Dc ,Rc ,dts,ts,dtl,tl) ; 

%  Evaluate  the  cost  function  separately  for  small  and  large  dt 
%  responses  and  sum  to  get  total  : 

r=[100];  %  indicates  desired  step  response 

jl  =  jcost(  xs,  r,  ys,  qmatrix,  rmatrix,  dts  ); 

j2  =  jcost(  xl,  r,  yl,  qmatrix,  rmatrix,  dtl  ); 

j  =  jl  +  j2  ; 

%  Evaluate  the  rise  time 


lens  =  length (  ts  )  ; 

t^rise  =  risetime(  xs(:.l),  ts,  .8  )  ; 

il  t_rise  >=  ts(  lens  ) 

t_rise  =  risetime(  xl(:,l),  tl,  ,8  )  ; 
end  ; 

%  Determine  infinity  norm  of  actuator  response 

act^norm  =  infnormC  xs(:,3)  )  ; 

%  Determine  the  Hang  Off  Error  in  state  1  (  settling  error  ) 

[  xrows,  xcols  ]  =  size(  xl  )  ; 

hang^off^error  =  xl(  xrows,  1  )  -  r(  1  )  ; 

if  isnan(  hang^off .error  ) 
hang.off .error  =  10  ; 
end  ; 

if  abs(  hang.off .error  )  >  10 
hang.off .error  =  10  ; 
end  ; 


%  Convert  cost  function  to  suitably  limited  fitness  function 
%  including  penality  for  exceeding  rise  time  specification 

if  isnan(  j  ) 
j  =  100  ; 
end  ; 

if  i  >  100 
j  =  100  ; 
end  ; 

if  t.rise  >  phi .rise 

cost .term  =  j  +  10.0  *  (  t.rise  -  phi.rise  )  ; 

else 

cost.term  =  j  ; 
end  ; 

%  Penalty  for  Settling  Error 

if  abs(  hang.off .error  )  >  hoe.thresh 

cost.term  =  cost.term  +  10.0  *  abs(  hang.off .error  ) 
end  ; 

%  Penalty  for  infinity  norm  of  actuator  response 

if  abs(  act .norm  )  >  act.thresh 
cost.term  =  cost.term  . . 

+  20.0  *  abs(  act.norm  -  act.thresh  ) 
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end  ; 

%  Limit  Cost  Term 

if  cost^term  >  100 
cost^term  =  100  ; 
elseif  cost^term  <  .6 
cost_term  =  .5  ; 

end  ; 

fitness  =  exp(  2  /  cost^term  )  -  1  ; 
pop^fitness  =  [  pop^fitness  ;  fitness  ]  ; 

%  Write  data  to  screen  if  not  in  continuous  plot  mode 
if  cont.plot  *"=  1 

k _ fitness  =  [  k  kg  fitness  ] 

end  ; 

end  7o  for  k  loop 


%  End  of  EVALP0P7 
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B.5  EVALPOPT 

This  MATLAB  function  evaluates  the  genetic  algorithm  fitness  function  for  the 
simple  single  parameter  example  used  in  program  GATEST.  It  is  only  used  with 
GATEST. 

A  program  listing  of  this  function  follows. 


y.  FUNCTION  EVALPOPT  -  Evaluate  Population  function  M  file 
%  -  assuming  one  parameter  per  string  :  kl 

%  -  special  function  for  GATEST  program 

% 

%  by  R.  Hull.  3/22/93 

function  [  pop_fitness,  pop_param  ]  =  .. 

evalpoptC  popsize,  kg_scale,  pop_chrom  ); 

pop_param  =  [  ]  ; 
pop_fitness  =  [  ]  ; 
kg  =  [  1.0  ]  : 

for  k  =  l:popsize  %  Population  loop 

kg(l)  =  conv2nm2(  pop_chrom(k,  1:16),  kg_scale(l) , . . 

kg_scale(2)  )  ; 

pop.param  =  [  pop.param  ;  kg  ]  ; 

%  Compute  the  special  fitness  function  of  one  parameter  : 
xl  =  kg(  1  )  ; 

fitness  =  -(  xl-2.5  )*  2  +  10.0  +  6*sin(  4*pi*xl  )  ; 

pop_fitness  =  [  pop_fitness  ;  fitness  ]  ; 

%  Write  data  to  screen  if  not  in  continuous  plot  mode 

if  cont.plot  ■=  1,  k _ fitness  =  [  k  kg  fitness  ],  end; 

end  %  for  k  loop 


y.  End  of  EVALPOPT 
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B.6  INFNORM 


This  MATLAB  function  computes  the  peak  response  value,  also  known  as  the  in¬ 
finity  norm,  of  the  given  input  signal. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  INFNORM 

%  This  function  computes  the  Infinity  Norm  of  a 
%  vector  function  input  y. 

% 

%  where  : 

%  y  =  input  response  vector 

%  norm  -  output  scaler  equal  to  the  infinity  norm 

%  of  y. 

% 

y.  written  in  MATLAB  by  R.  Hull  2/18/93 


function  norm  =  inf norm (  y  ) 
abs_y  =  abs(  y  )  :  % 

norm  =  max(  abs_y  )  ;  % 

%  End  of  function  infnorm 


a  vector 
a  scaler 
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B.7  JCOST 


This  MATLAB  function  computes  the  weighted  linear  quadratic  performance  index, 
given  the  time  history  of  the  system  state  and  control  variables. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  JCOST 

%  This  function  computes  the  Linear  Quadratic  Requlator 
%  cost  function  based  on  matrices  containing  the 
%  time  history  of  the  input  and  states  of  the  system 
%  where  : 

%  1  =  number  of  time  samples 

%  n  =  number  of  states 

%  m  =  number  of  control  inputs 

%  X  =  history  of  the  states  (  1  x  n  ) 

%  r  =  desired  step  response  of  states  (  1  x  n  ) 

%  u  =  history  of  the  controls  (  1  x  m  ) 

%  q  =  state  weighting  matrix  (  n  x  n  ) 

%  R  =  control  weighting  matrix  (  m  x  m  ) 

%  dt  =  time  sample  interval  (  constant  ) 

%  written  in  MATLAB  by  R.  Hull  10/14/92 

function  [j]  =  jcost(x,r ,u,Q,R,dt) 

j  =  0  : 
jdotp  =  0  ; 

for  i  =  1 : length (x( : ,1))  ; 
xerr  =  r  -  x(i , : )  ; 

jdot  =  xerr  *  Q  *  xerr.’  +  u(i,:)  *  R  *  u(i,;),’  ; 
j=i+.5*dt*(  jdot  +  jdotp  )  ; 
jdotp  =  jdot  ; 
end; 

%  End  of  Function  JCOST 
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B.8  RISETIME 


This  MATLAB  function  computes  risetime,  given  a  response  history,  as  measured 
from  0  to  80  %  of  the  commanded  value. 

A  program  listing  of  this  function  follows. 


% 

% 

% 

7. 

7. 

7. 

7. 

7. 

% 

7. 

7. 

7. 

7. 

7. 

7. 

7. 

7. 


FUNCTION  RISETIME 

This  function  computes  the  Rise  Time  of  a  step 
response  function  to  the  value  input  as  y_value. 

If  the  step  response  never  reaches  the  desired 
y_value,  the  rtime  returned  will  be  the  final 
value  in  the  time  vector. 

where  : 

y  =  step  response  vector 
t  =  time  vector 

y_value  =  desired  value  of  y  at  which  response  time 
is  measured 

r_time  =  earliest  time  for  which  step  response  in  y 
always  exceeds  y_value.  Note  that  no 
interpolation  is  performed. 

written  in  MATLAB  by  R.  Hull  1/25/93 


function  rtime  =  risetime (  y,  t,  y_value  )  ; 

y_index  =  length (  y  )  ; 

for  i  =  1  :  lengthC  y  ) 
if  y(i)  >=  y_value 
y_index  =  i  ; 
break  ; 
end  ; 
end  ; 

len_t  =  lengthC  t  )  ; 

if  y_index  <  len_t 

rtime  =  t(  y_index  )  ; 
else 

rtime  =  t(  len_t  )  ; 

end  ; 

%  End  of  function  RISETIME 
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B.9  SELECT 

This  MATLAB  function  selects  two  parents  from  the  available  population  pool 
for  mating  as  part  of  the  genetic  algorithm.  The  probability  of  selection  is  based 
on  a  uniformly  distributed  random  function,  but  is  weighted  according  to  each 
individual’s  fitness  relative  to  the  total  population  fitness.  It  ensures  that  the  two 
parents  selected  are  different  individuals. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  SELECT  -  Select  Parents  for  Genetic  Algorithm 

% 

i  by  R.  Hull  1/12/93 

% 

%  Note  :  this  selection  procedure  insures  that 
%  parent_l  and  parent_2  are  different  individuals. 

function  [  parent_l,  parent_2  ]  =  select (  popsize,  sum^fitness,  .. 

pop_fitness  )  ; 

%  Determine  random  roulette  wheel  pointer 
pointer_sum  =  rand  *  sum^fitness  ; 
part_sum  =  0  ; 
for  parental  =  1  :  popsize 

part_sum  =  part^sum  +  pop_fitness(  parental  )  ; 
if  part_sum  >=  pointer^sum,  break,  end  ;  %  break  loop 

end  ;  %  (  for  parental  loop  ) 

%  Remove  parental  from  choices  and  repeat  to  find  parent_2 
sum.fitness  =  sum^fitness  -  pop_fitness(  parent.l  )  ; 
pop_fitness(  parental  )  =  0  ; 

%  Determine  second  random  roulette  wheel  pointer 
pointer^sum  =  rand  *  sum^fitness  ; 
part_sum  =0  ; 
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for  pareiit_2  =  1  :  popsize 

part_sum  ”  part^sum  +  pop_f itnessC  parent_2  )  ; 
if  part^sum  >=  pointer^sum,  break,  end  ;  %  break  loop 

end  ;  %  (  for  parent_2  loop  ) 

%  End  of  function  SELECT 
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B.IO  SELXSITE 


This  MATLAB  function  selects  the  crossover  site  for  the  crossover  function  of  the 
genetic  algorithm.  It  implements  an  experimental  crossover  function,  in  which  the 
random  distribution  of  the  site  selection  is  weighted  in  favor  of  high  order  bits  in 
early  generations,  gradually  shifting  to  favor  low  order  bits  in  later  generations. 

A  program  listing  of  this  function  follows. 


%  FUNCTION  SELXSITE  -  Select  Genetic  Algorithm  Crossover  Site 

% 

%  by  R.  Hull  2/15/93 

% 

%  This  function  selects  a  crossover  site  for  the 

%  Genetic  Algorithm  that  is  stochastically  weighted 

%  depending  upon  generation. 

%  It  uses  a  gaussian  distribution  that  is  "folded  over" 

%  at  the  ends  of  the  parameter  bit  length. 


function  xsite  =  selxsiteC  paramlength,  generation  )  ; 

%  Switch  to  Normal  Distribution 

rand ( ’normal’ )  ; 

sigma  =  paramlength  /  3  ; 

mean  =  generation  ; 

if  mean  <  1 
mean  =  1  ; 

else if  mean  >  paramlength 
mean  =  paramlength  ; 
end  ; 

gauss.max  =  paramlength  -  1  ; 

7*  Generate  Gaussian  number  and  fold  over  negative  tail 

gauss^num  =  abs(  rand  ♦  sigma  +  mean  )  ; 

%  Fold  over  positive  side  tail 

if  gauss.num  >  (  2  *  gauss^max  ) 
gauss.num  =  2  *  gauss_max  ; 
end  ; 


±f  gauss_num  >  gauss_roax 

gauss^nura  =  2  *  gauss^max  ~  gauss^num  ; 
end  ; 

%  Convert  to  integer  crossover  site  between  1  and  paramlength 
xsite  =  iix(  gauss.num  )  +  1  ; 

%  Switch  Back  to  Uniform  Distribution 
rand ( ’uniform* )  ; 

%  End  of  function  SELXSITE 
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B.ll  STEP2S 

This  MATLAB  function  computes  the  step  response  to  a  LTIC  system  using  two 
different  integration  time  intervals.  It  first  uses  a  small  dt  (dts)  to  provide  accurate 
integration  during  the  initial  period  of  high  dynamics.  It  then  switches  to  a  larger 
dt  (dtl)  to  compute  the  steady  state  step  response  during  a  period  of  lower  dynamic 
rates.  This  enables  an  accurate  response  to  be  generated,  yet  conserves  computation 
time  and  resources. 

A  program  listing  of  this  function  follows. 


y.  Function  STEP2S 

%  Computes  the  step  response  of  continuous-time  linear 

%  systems,  using  split  dt's. 

% 

%  This  function  calculates  the  response  of  the  system: 

y.  . 

y,  X  =  Ax  +  Bu 

y.  y  =  Cx  +  Du 

y. 

%  to  a  step  input  vector  R  (  m  x  1  ) .  using  dts  over 

%  the  time  history  in  ts,  and  dtl  over  the  time  history 

%  in  tl . 

% 

%  ys  is  the  response  over  time  history  ts 

%  yl  is  the  response  over  time  history  tl 

%  xs  is  the  state  history  over  time  history  ts 

%  xl  is  the  state  history  over  time  history  tl 

y. 

%  This  allows  an  accurate  step  response  to  be  determined  over 

%  the  initial  period  of  high  actuator  dynamics  using  a  small 

y,  dt  (dts) .  then  a  steady  state  response  to  be  generated  using 

y,  a  large  dt  (dtl)  ,  without  using  too  much  computer  time . 

%  The  final  states  of  the  small  dt  response  become  the  initial 

y,  states  for  the  large  dt  response. 

i  by  R.  Hull  3/1/93 

function  [ys,yl,xs,xl]  =  step28(a,b,c.d,r,dts,ts,dtl,tl)  ; 

%  Convert  Continuous  System  to  Discrete  Time  Systems  : 

[aas.bbs]  =  c2d(a,b,dts) ; 

[aal.bbl]  =  c2d(a.b.dtl) ; 


%  Build  Small  Step  Input  Matrix  from  input  vector 
lens  -  length(ts) ; 
u  =  ones(  lens,  1  )  *  r. *  ; 

%  Compute  linear  time  invariant  state  response  : 
xs  =  ltitr(aas,bbs,u) ; 

%  Compute  output  vector  : 
ys  =  xs  *  c . •+  u  *  d. * ; 

%  Build  Large  Step  Input  Matrix  from  input  vector 
lenl  =  length (tl) ; 
u  =  ones(  lenl,  1  )  *  r.*  ; 

%  Initialize  states  to  final  values  from  small  dt  response 
xsf  =  xs(  lens,  :  )  ; 

%  Compute  linear  time  invariant  state  response  : 
xl  =  ltitr(aal,bbl,u,xsf ) ; 

%  Compute  output  vector  : 
yl=xl*c.*+u*d.'; 

7o  End  of  function  STEP2S 


